Sagar S Vyavahare

Aleksandra Mieczkowska

Peter R Flatt

Daniel Chappard

Nigel Irwin

PhD Guillaume Mabilleau
email: guillaume.mabilleau@univ-angers.fr

GIP analogues augment bone strength by modulating bone composition in diet-induced obesity in mice

Keywords: GIP-Tag, bone fragility, bone composition, type 2 diabetes

D-Ala²)-GIP and (D-Ala²)-GIP-Tag had similar biological effects in vitro • (D-Ala²)-GIP and (D-Ala²)-GIP-Tag reduced non-

INTRODUCTION

The burden of diabetes is progressing worldwide at an alarming rate. It is projected that by 2045 there will be an estimated number of 629 million individuals with diabetes [START_REF] Federation | IDF Diabetes Atlas[END_REF]. Diabetes is associated with comorbidities including micro-and macrovascular complications, peripheral neuropathy and bone fractures. Indeed, diabetes is associated with a higher risk of bone fractures at the hip and at the spine, despite higher or normal bone mineral densities at these two sites [START_REF] Schwartz | Study of Osteoporotic Features Research, Older women with diabetes have an increased risk of fracture: a prospective study[END_REF][START_REF] Vestergaard | Discrepancies in bone mineral density and fracture risk in patients with type 1 and type 2 diabetes--a meta-analysis[END_REF][START_REF] Bonds | Risk of fracture in women with type 2 diabetes: the Women's Health Initiative Observational Study[END_REF][START_REF] Napoli | Osteoporotic Fractures in Men Study Research, Vertebral Fracture Risk in Diabetic Elderly Men: The MrOS Study[END_REF].

The resistance of bone to fracture depends on bone mass, i.e., the quantity of bone, but also on the organization of the bone matrix at the structural and tissue levels. Alterations of both cortical and trabecular bone microarchitectures have been linked to decreased bone strength. Bone composition is also an important contributor of bone strength. The bone matrix is comprised of an organic layer, composed mainly of type I collagen, proteoglycans and non-collagenous protein [START_REF] Boskey | The composition of bone[END_REF]. This organic layer serves as a template for the nucleation and growth of poorly crystalline hydroxyapatite and the subsequent mineralization of bone matrix. Both the organic and mineral layers are important for bone strength, and alterations in either of these components leads to changes in bone strength.

Glucose-dependent insulinotropic polypeptide (GIP) is produced and secreted by enteroendocrine Kcells, located mostly in the duodenum [START_REF] Habib | Overlap of endocrine hormone expression in the mouse intestine revealed by transcriptional profiling and flow cytometry[END_REF]. Whilst in the circulation, GIP can access several target organs such as the pancreas where it stimulates insulin secretion from beta pancreatic cells [START_REF] Baggio | Biology of incretins: GLP-1 and GIP[END_REF]. Physiologically GIP is one of the two incretin hormones responsible for potentiation of insulin secretion following feeding [START_REF] Nauck | Incretin hormones: Their role in health and disease[END_REF]. We previously reported beneficial effects of several gut hormones on bone strength in diabetes [START_REF] Mansur | A new stable GIP-Oxyntomodulin hybrid peptide improved bone strength both at the organ and tissue levels in genetically-inherited type 2 diabetes mellitus[END_REF][START_REF] Mansur | Sitagliptin Alters Bone Composition in High-Fat-Fed Mice[END_REF][START_REF] Mansur | The GLP-1 Receptor Agonist Exenatide Ameliorates Bone Composition and Tissue Material Properties in High Fat Fed Diabetic Mice[END_REF][START_REF] Pereira | Exenatide Improves Bone Quality in a Murine Model of Genetically Inherited Type 2 Diabetes Mellitus[END_REF]. However, the possible effects of GIP analogues have not been reported to date. Due to its role as the major physiological incretin hormone, GIP effects on bone are interesting to test in diabetes.

Deletion of GIP receptor (GIPr) in mice leads to a skeletal phenotype with deterioration of cortical bone microarchitecture, alterations of both mineral and collagen component of the bone matrix and ultimately lower bone strength [START_REF] Gaudin-Audrain | Glucose-dependent insulinotropic polypeptide receptor deficiency leads to modifications of trabecular bone volume and quality in mice[END_REF][START_REF] Mieczkowska | Glucose-dependent insulinotropic polypeptide (GIP) receptor deletion leads to reduced bone strength and quality[END_REF][START_REF] Tsukiyama | Gastric inhibitory polypeptide as an endogenous factor promoting new bone formation after food ingestion[END_REF][START_REF] Xie | Glucose-dependent insulinotropic polypeptide receptor knockout mice have altered bone turnover[END_REF]. Administration of the long-acting, N-AcGIP analogue in rats, led to ameliorations of bone strength due to improvement in bone composition [START_REF] Mabilleau | Beneficial effects of a N-terminally modified GIP agonist on tissue-level bone material properties[END_REF]. The GIPr is expressed by bone cells (osteoblasts and osteoclasts) and culture of GIP analogues with bone cells enhances osteoblast function and reduces osteoclast formation and activation [START_REF] Mieczkowska | Glucose-dependent insulinotropic polypeptide (GIP) directly affects collagen fibril diameter and collagen cross-linking in osteoblast cultures[END_REF][START_REF] Mabilleau | Glucosedependent insulinotropic polypeptide (GIP) dose-dependently reduces osteoclast differentiation and resorption[END_REF]. In humans, GIP also exerts anti-resorptive effects [START_REF] Nissen | Glucose-dependent insulinotropic polypeptide inhibits bone resorption in humans[END_REF][START_REF] Christensen | Glucose-Dependent Insulinotropic Polypeptide (GIP) Inhibits Bone Resorption Independently of Insulin and Glycemia[END_REF] and a polymorphism of the GIPr has been linked to higher risk of developing non vertebral bone fractures [START_REF] Torekov | A functional amino acid substitution in the glucose-dependent insulinotropic polypeptide receptor (GIPR) gene is associated with lower bone mineral density and increased fracture risk[END_REF].

GIPr are found in extra-skeletal tissues, and it is therefore unclear whether the beneficial effects on bone are due to direct activation of bone cells or by indirect action on extra-skeletal organs that are responsible for better bone strength. To overcome this problem, we developed a novel GIP analogue that exhibits specific bone tropism [START_REF] Mabilleau | Efficacy of targeting bonespecific GIP receptor in ovariectomy-induced bone loss[END_REF]. This analogue, called GIP-Tag, has N-terminal enzymatic protection by substitution of L-Ala 2 with D-Ala 2 , together with six aspartic acids at its C-terminal extremity that confers affinity for bone mineral. We have previously shown that GIP-Tag exhibits comparable binding activity at the GIPr and activates similar intracellular pathways in bone cells as the parent GIP peptide [START_REF] Mabilleau | Efficacy of targeting bonespecific GIP receptor in ovariectomy-induced bone loss[END_REF].

Recent works highlighted the dramatic bone fragility resulting from high fat diet on bone strength and quality [START_REF] Mansur | Sitagliptin Alters Bone Composition in High-Fat-Fed Mice[END_REF][START_REF] Mansur | The GLP-1 Receptor Agonist Exenatide Ameliorates Bone Composition and Tissue Material Properties in High Fat Fed Diabetic Mice[END_REF] and as such, the aim of the present study was to evidence whether GIP-Tag was capable of enhancing osteoblast mineralization and collagen post-processing in vitro, and to assess whether GIP-Tag was capable of ameliorating bone strength in a diet-induced obesity rodent model.

MATERIAL AND METHODS

Peptide characterization

(D-Ala 2 )-GIP and (D-Ala 2 )-GIP-Tag were synthesized and purchased from EZ Biolabs Ltd (Carmel, IN, USA) at >95% purity. Peptides were further characterized in-house by reverse phase high performance liquid chromatography after dissolution in water using a Kinetex C-18 HPLC analytical column (150 x 4.6 mm, Phenomenex, Cheshire, UK). Equilibration of the column was performed with 0.1% (v/v) TFA/water at a flow rate of 1 ml/min. Linear acetonitrile gradients (0-70% in 40 mins) were used for elution and absorbance measured at 214 nm. The molecular mass of peptides was confirmed using a Voyager-DE biospectrometry workstation (PerSeptive Biosystems, Farmingham, MA, USA). Dipeptidyl peptidase-4 (DPP-4) mediated degradation of the peptides was assessed as described by Gault et al [START_REF] Gault | Enhanced cAMP generation and insulin-releasing potency of two novel Tyr1-modified enzyme-resistant forms of glucose-dependent insulinotropic polypeptide is associated with significant antihyperglycaemic activity in spontaneous obesity-diabetes[END_REF]. Briefly, 30 µg peptides were incubated with 5 mU porcine DPP-4 in 50 mM triethanolamine hydrochloride (TEA, pH 7.8) at 37°C for 0, 2, 8, 12 and 24 h. The enzymatic reaction was then stopped using 50 µl TFA/water (10% v/v). Degradation profiles were determined using HPLC and degradation products were analyzed by mass spectrometry.

Cell culture

MC3T3-E1 cells are phenotypically normal osteoblasts derived from foetal mouse calvaria and differentiate in culture resulting ultimately to the formation of a bone-like extracellular matrix [START_REF] Franceschi | Effects of ascorbic acid on collagen matrix formation and osteoblast differentiation in murine MC3T3-E1 cells[END_REF].

MC3T3-E1 cells were purchased from American Type Culture Collection (ATCC, Teddington, UK). Cells were grown and expanded at a ratio 1:4 in propagation medium containing αMEM supplemented with 10% FBS, 100 U/mL penicillin, and 100 μg/mL streptomycin in a humidified atmosphere enriched with 5% CO2 at 37°C. For differentiation studies, cells were detached with trypsin-EDTA and plated at a density of 1.5 × 10 4 cells/cm 2 and grown to confluence in propagation medium. At confluence, the propagation medium was replaced by differentiation medium containing αMEM supplemented with 10% FBS, 100 U/mL penicillin, 100 μg/mL streptomycin, 3.5 g/L D-glucose (final glucose concentration 4.5 g/L), 50 µg/ml ascorbic acid, 2.5 mM Na2HPO4 and 100 pM peptide. The differentiation medium was replenished every two days. At the end of the 13-day culture period, cells were either fixed with 3.7% (w/v) formaldehyde in PBS for 10 minutes (mineralization assay) or absolute ethanol (collagen maturity assay). For the mineralization assay, the cell layer was then rinsed with PBS and stained with 1% (w/v) alizarin red S for 5 minutes. The plates were then rinsed three times with 50% ethanol and left to air dry.

An image of the full well was made by light microscopy by inverting the plate on the microscope stage.

The extent of mineralization was assessed by image analysis with ImageJ 1.52p (National Institutes of Health, Bethesda, MA, USA). A lower grey level indicates a higher mineralization.

Lysyl oxidase (LOX) activity

LOX activity was assessed in the cell culture supernatants and in the cell layers after 13 days of culture. This period of time was chosen as it corresponds to the highest expression of LOX [START_REF] Hong | A role for lysyl oxidase regulation in the control of normal collagen deposition in differentiating osteoblast cultures[END_REF] using a fluorometric assay developed by Palamakumbura et al. [START_REF] Palamakumbura | A fluorometric assay for detection of lysyl oxidase enzyme activity in biological samples[END_REF]. At day 12, the differentiation medium was replaced by phenol red-free DMEM containing 0.5% bovine serum albumin, 50 µg/ml ascorbic acid, 2.5 mM Na2HPO4 and 100 pM peptide. After 24 hrs, cell culture supernatants were collected, centrifuged at 3000 rpm for 30 min at 4°C, aliquoted and stored at -80°C until use. The cell layer was extracted in urea buffer as reported previously [START_REF] Bedell-Hogan | Oxidation, cross-linking, and insolubilization of recombinant tropoelastin by purified lysyl oxidase[END_REF]. Samples (supernatant, cell layer extract) were then incubated at 37°C for 30 minutes in the presence of 1.2 M urea, 50 mM sodium borate (pH 8.2), 1.3 nmol H2O2, 1 UI/ml horseradish peroxidase, 10 mM diaminopentane, 10 µM Ampliflu TM red and ± 0.5 mM betaaminopropionitrile (BAPN) in opaque 96-well plates. At the end of the incubation period, plates were placed on ice and fluorescence was read using a M2 microplate reader (Molecular devices, St Gregoire, France) with excitation and emission wavelengths set up at 563 nm and 587 nm, respectively. As BAPN is a specific LOX inhibitor, the residual amine oxidase activity evidenced in the presence of BAPN was subtracted from the activity observed in the absence of BAPN. LOX activity was reported as the fold change versus CTRL of the pooled supernatant and cell layer fractions.

Osteoclast generation from human peripheral blood mononuclear cells

Peripheral blood mononuclear cells (PBMCs) were isolated from the buffy coat of healthy individuals (Etablissement français du sang, Angers, France) as reported previously [START_REF] Mabilleau | Role of the A20-TRAF6 axis in lipopolysaccharidemediated osteoclastogenesis[END_REF]. In order to generate osteoclasts, PBMCs cells were plated at a concentration of 1 x 10 6 cells/cm² in 24-well plate or on top of dentine slices in alpha-MEM supplemented with 10% fetal calf serum, 2 mM L-glutamine, 100 U/ml penicillin, 100 µg/ml streptomycin, 50 ng/ml sRANKL (Bio-techne) and 10 nM of GIP analogues. This concentration of GIP was used based on previous published data [START_REF] Mabilleau | Glucosedependent insulinotropic polypeptide (GIP) dose-dependently reduces osteoclast differentiation and resorption[END_REF]. All factors were replenished every 2-3 days. After 14 days of culture, the expression of tartrate resistant acid phosphatase (TRAP) was examined histochemically in 24-well plate. Briefly, cells were rinsed promptly in PBS buffer, fixed with formalin (10% in PBS buffer) for 10 minutes and rinsed in distilled water. TRAP was demonstrated histochemically by a simultaneous coupling reaction using Naphtol AS-BI-phosphate as substrate and Fast violet B as the diazonium salt. Cells were then incubated for 90 minutes at 37°C in the dark, rinsed three times in distilled water and the residual activity was inhibited by 4% NaF for 30 minutes. Cells were then rinsed in distilled water, counterstained with DAPI for 20 minutes and allowed to dry. TRAP positive cells, with more than three nuclei, were identified as osteoclasts. Furthermore, after 21 days in culture, the dentine slices were removed from the culture wells, placed in NH4OH (1N) for 30 minutes and sonicated for 5 minutes to remove any adherent cells. After rinsing in distilled water, the dentine slices were stained with 0.5% (v/v) toluidine blue (pH 5.0) prior to examination by light microscopy. The extent of surface erosion were determined using imageJ.

Animals

All procedures were conducted in accordance with UK Home Office Regulations (UK Animals Scientific Procedures Act 1986) and in compliance with the ARRIVE guidelines [START_REF] Kilkenny | Improving bioscience research reporting: the ARRIVE guidelines for reporting animal research[END_REF]. Male NIH Swiss mice (NIH/OlaHsd) were obtained from Envigo Ltd (Blackthorn, UK) at 8 weeks of age. Animals were individually housed in an air-conditioned room at 22 ± 2°C with a 12-hour light/dark cycle, and had free access to tap water and high fat diet (45% fat, 20% protein, 35% carbohydrate; Special Diet Service, Essex, UK) ad libitum for 90 days prior to the start of the study. High fat fed (HFF) mice (n=6) were divided in three groups and received once daily (09:30 h) intraperitoneal (i.p.) injections of: (1) saline (0.9% NaCl -HFD+saline), (2) (D-Ala²)-GIP (25 nmol/kg bw) and (3) (D-Ala²)-GIP-Tag (25 nmol/kg bw), with both peptides delivered in saline vehicle. This dosing regimen was chosen based on our previous experience with related gut-derived peptides in high fat fed mice [START_REF] Mansur | The GLP-1 Receptor Agonist Exenatide Ameliorates Bone Composition and Tissue Material Properties in High Fat Fed Diabetic Mice[END_REF][START_REF] Irwin | Metabolic effects of sustained activation of the GLP-1 receptor alone and in combination with background GIP receptor antagonism in high fat-fed mice[END_REF][START_REF] Martin | A novel acylated form of (d-Ala(2))GIP with improved antidiabetic potential, lacking effect on body fat stores[END_REF]. Pharmacological intervention was continued for 42 consecutive days. No adverse effects were observed following peptide treatment. Body weight and blood glucose concentrations were measured at regular intervals, blood was withdrawn at approximately 10:00 h on appropriate days. An i.p. glucose tolerance test (18 mmol/kg bw) was performed following an 18 hour fast on day 43. All blood samples were collected from the tail vein of animals and glucose measured using a portable Bayer Ascencia Counter blood glucose meter (Bayer Healthcare, Newbury, Berkshire, UK). At necropsy, femurs and tibias were cleaned of soft tissues and stored in 70% ethanol at 4°C until used.

MicroCT

MicroCT analyses were performed on tibias with a Skyscan 1272 microtomograph (Bruker MicroCT, Kontich, Belgium) operated at 70 kV, 140 µA, 1000-ms integration time. The isotropic pixel size was fixed at 4 µm, the rotation step at 0.25° and exposure was performed with a 0.5-mm aluminum filter.

Each 3D reconstruction image dataset was binarized using global thresholding. Trabecular volume of interest (VOI) was located 0.5-mm below the growth plate and extended on 2-mm. Cortical volume of interest extended on 1-mm centered at the midshaft tibia. All histomorphometrical parameters were determined using the CTan software (version 1.17.7.2, Bruker microCT) according to guidelines and nomenclature proposed by the American Society for Bone and Mineral Research [START_REF] Bouxsein | Guidelines for assessment of bone microstructure in rodents using micro-computed tomography[END_REF].

Bone strength assessment

Three-point bending experiments were performed on femurs after rehydrating bones at 4°C for 24 h. Femurs were loaded to failure in 3-point bending at 1 mm/min using an Instron 5942 (Instron, Elancourt, France). The lower span length was set to 10 mm. Femurs were oriented so the anterior quadrant was facing down and subjected to tensile loads. For each bone, load and displacement were digitally recorded at a sampling rate of 100 Hz and measured using a 500 N load cell (Instron). The loaddisplacement curve was computed with the Bluehill 3 software (Instron) and the maximum load, yield load (0.2% offset method), stiffness, post-yield displacement and work to failure were computerized.

After three-point bending experiments, distal parts of femur were embedded undecalcified in polymethylmethacrylate (pMMA) at 4°C.

Fourier-transform infrared microspectroscopy/imaging (FTIRM/FTIRI)

For collagen cross-link analysis (collagen maturity), osteoblast cultures were fixed in absolute ethanol as indicated above, scrapped off the culture dish and transferred onto BaF2 windows where they were air-dried. Integrity of the collagen extracellular matrix was verified by comparing the obtained FTIRM spectrum with those of commercial collagen. Spectral analyses were performed using a Bruker Vertex 70 spectrometer (Bruker optics, Ettlingen, Germany) interfaced with a Bruker Hyperion 3000 infrared microscope equipped with a standard single element Mercury Cadmium Telluride (MCT) detector.

Infrared spectra were recorded at a resolution of 4 cm -1 , with an average of 32 scans in transmission mode. Background spectral images were collected under identical conditions from the same BaF2 windows at the beginning and end of each experiment to ensure instrument stability. For FTIRM analysis, 20 spectra were acquired for each condition. Each spectrum was corrected for Mie scattering with the RMieS-EMSC_v5 algorithm (kind gift of Prof Peter Gardner, University of Manchester, UK). Spectra were vector normalized for the amide I band, second derivative spectroscopy was applied to find the position of underlying peaks in the amide I band and curve fitting was performed with a routine script in Matlab (The Mathworks, Natick, USA). Collagen maturity was determined as the relative ratio of pyridinium trivalent (Pyr, mature collagen) to dehydrodihydroxylysinonorleucine divalent (deH-DHLNL, new collagen) collagen cross-links using their respective subbands located at 1660 cm -1 and 1690 cm -1 of the amide I peak [START_REF] Paschalis | Fourier Transform Infrared Spectroscopic Characterization of Mineralizing Type I Collagen Enzymatic Trivalent Cross-Links[END_REF][START_REF] Paschalis | Spectroscopic characterization of collagen cross-links in bone[END_REF].

For analysis of in vivo samples, cross-sectional sections (4 µm) of the midshaft femur were sandwiched between barium fluoride optical windows. FTIRI was performed with the same spectrometer and microscope but with a focal plane array detector (64 x 64 pixels) covering a field of view of 180 x 180 µm. Nine consecutive field-of-views were stitched together to allow sufficient bone to be analyzed.

Sections were scanned with a spectral resolution of 8 cm -1 (spectral region 900-2000 cm -1 ). Each spectrum was corrected for Mie scattering. Prior to pMMA subtraction, the area of the 1715 -1745 cm - cm -1 were computed. After pMMA subtraction, spectrum were vector normalized for the v1,v3 phosphate band analyzed with a routine script in Matlab (The Mathworks, Natick, USA) as previously reported [START_REF] Aguado | Hypodynamia Alters Bone Quality and Trabecular Microarchitecture[END_REF].

The evaluated infrared spectral parameters were: (1) Tissue water content, calculated as the area ratio of the pMMA peak (1715-1745 cm -1 ) to amide I (1590-1710 cm -1 ), both computed before pMMA subtraction. This ratio represents a surrogate to water content in the bone matrix [START_REF] Paschalis | Aging Versus Postmenopausal Osteoporosis: Bone Composition and Maturation Kinetics at Actively-Forming Trabecular Surfaces of Female Subjects Aged 1 to 84 Years[END_REF]. (2) mineral/matrix ratio (area of υ1, υ3 phosphate/area amide1); (3) acid phosphate content (intensity ratio 1127cm -1 /1096 cm -1 ); (4) mineral maturity (intensity ratio 1030 cm -1 /1020 cm -1 ), reflecting crystal size and perfection; [START_REF] Napoli | Osteoporotic Fractures in Men Study Research, Vertebral Fracture Risk in Diabetic Elderly Men: The MrOS Study[END_REF] crystal size index (intensity ratio 1075 cm -1 /1055 cm -1 ), reflecting crystal size in 002, 211, 200 and 202 directions [START_REF] Gadaleta | Fourier transform infrared spectroscopy of the solution-mediated conversion of amorphous calcium phosphate to hydroxyapatite: new correlations between X-ray diffraction and infrared data[END_REF]; (6) carbonate/phosphate ratio (intensity υ3 carbonate located at ~1415 cm -1 /1030 cm -1 ) was computed after subtracting the organic matrix spectrum [START_REF] Ou-Yang | Infrared microscopic imaging of bone: spatial distribution of CO3(2-)[END_REF]; and (7) collagen maturity (intensity ratio 1660 cm -1 /1690 cm -1 ). Histogram distribution for each compositional parameter was fitted with a Gaussian model and considered normally distributed if the R² coefficient was >0.95. In the present study, none of the histograms deviated from normal distribution. For each of the compositional parameters, the mean and full width at half maximum of the pixel distribution (excluding the zero background values) were computed and represented as mean and heterogeneity.

Statistical analysis

All data were analyzed using Prism 6.0 (GraphPad Software Inc, La Jolla, CA). Kruskal-Wallis with Dunn's multiple comparisons tests were employed to test for significance when otherwise stated. Twoway ANOVA with Tukey multiple comparison tests were used to assess significant differences in body weight and blood glucose. Multiple regression analyses were performed with Systat 13 (Systat software, San Jose, CA, USA). Differences at p equal to or less than 0.05 were considered significant. Data are presented as median ± interquartile range.

RESULTS

(D-Ala 2 )-GIP-Tag and (D-Ala 2 )-GIP acted similarly on bone cells in vitro

Initially, as expected, we confirmed that both (D-Ala²)-GIP and (D-Ala²)-GIP-Tag were resistant to DPP-4 mediated degradation (data not shown) as both peptides were still intact after a 8-h incubation with DPP-4. We previously observed that (D-Ala²)-GIP and (D-Ala²)-GIP-Tag exhibited similar bone related responses in the presence of 1 g/L glucose [START_REF] Mabilleau | Efficacy of targeting bonespecific GIP receptor in ovariectomy-induced bone loss[END_REF]. However, little was known about the response in the presence of elevated glucose concentration. We assessed the extent of mineralization in the presence of 4.5 g/L glucose. As seen in Figure 1A, (D-Ala²)-GIP and (D-Ala²)-GIP-Tag significantly increased mineral deposition in the cultures. Furthermore, (D-Ala²)-GIP and (D-Ala²)-GIP-Tag significantly increased LOX activity, responsible for enzymatic collagen crosslinking (Fig. 1B). Indeed, collagen maturity, expressed as the ratio of mature/immature collagen crosslinks, was significantly augmented by 42% (p=0.026) and 43% (p=0.021) in (D-Ala²)-GIP-and (D-Ala²)-GIP-Tag-treated cultures, respectively. We next examined whether (D-Ala²)-GIP and (D-Ala²)-GIP-Tag were capable of reducing human osteoclast formation and activity in vitro. Human peripheral mononuclear cells were differentiated in the presence of (D-Ala²)-GIP or (D-Ala²)-GIP-Tag. As shown in Figures 1C-D, (D-Ala²)-GIP was capable of reducing osteoclast formation by 24% (p=0.022) and activity by 26% (p=0.014). Similarly, (D-Ala²)-GIP-Tag reduced osteoclast formation by 27% (p=0.027) and osteoclast-mediated bone resorption by 32% (p=0.018). No significant differences between the two pharmacological interventions were evident in any of the in vitro parameters studied.

(D-Ala²)-GIP and (D-Ala²)-GIP-Tag improved blood glucose and bone strength

As (D-Ala²)-GIP and (D-Ala²)-GIP-Tag were both active in vitro, we decided to administer these two peptides to high fat fed (HFF) animals. Previously, we reported that HFF animals presented with alterations of bone mechanical resistance, as well as deterioration not only in trabecular and cortical bone microarchitectures, but also in the mineral and organic components of bone matrix [START_REF] Mansur | Sitagliptin Alters Bone Composition in High-Fat-Fed Mice[END_REF][START_REF] Mansur | The GLP-1 Receptor Agonist Exenatide Ameliorates Bone Composition and Tissue Material Properties in High Fat Fed Diabetic Mice[END_REF]. HFF mice were characterized by increased body weight and elevated blood glucose. We first investigated therefore whether (D-Ala²)-GIP and (D-Ala²)-GIP-Tag affected these parameters. Body weight was not significantly modified during the 42 days administration of (D-Ala²)-GIP or (D-Ala²)-GIP-tag (Figure 2A). However, non-fasting blood glucose levels were significantly reduced after 42 days of treatment with (D-Ala²)-GIP and (D-Ala²)-GIP-Tag, respectively (Figure 2B). Interestingly, glucose levels following a glucose challenge on day 42 were reduced in (D-Ala²)-GIP (p<0.001) and (D-Ala²)-GIP-Tag (p=0.001) mice at 105 min post-injection, with glucose-lowering effects more prominent in (D-Ala²)-GIP treated mice (Figure 2 C-D). Corresponding glucose-induced insulin concentrations were not significantly different between all groups of HFF mice (data now shown). We also examined whether bone mechanical response was modified by peptide administration in HFF mice. Ultimate load was significantly increased by 13% (p=0.024) and 10% (p=0.041) with (D-Ala²)-GIP and (D-Ala²)-GIP-Tag, respectively (Figure 3). Work-to-fracture, representing the energy required to break the bone, was significantly augmented by 30% (p=0.004) and 20% (p=0.037) with (D-Ala²)-GIP and (D-Ala²)-GIP-Tag, respectively. On the other hand, stiffness and post-yield displacement were not significantly different between the three groups of animals.

(D-Ala²)-GIP and (D-Ala²)-GIP-Tag had no effects on bone microarchitecture

As mechanical response depends on structural and compositional properties, we assessed trabecular and cortical bone microarchitectures. Figure 4A represents 3D-models of tibia proximal metaphysis in saline-, (D-Ala²)-GIP-and (D-Ala²)-GIP-Tag-treated animals. No obvious differences could be detected between the three groups. Quantification of histomorphometrical parameters did not evidence any significant modifications of either trabecular or cortical bone microarchitecture (Figures 4B-C).

(D-Ala²)-GIP and (D-Ala²)-GIP-Tag effects on bone composition

We next examined whether pharmacological intervention with (D-Ala²)-GIP and (D-Ala²)-GIP-Tag modified the composition of the bone matrix. Each compositional parameters is described by its average quantity (mean) and heterogeneity (width). As reported in Figure 5, (D-Ala²)-GIP significantly augmented XLR-mean, XLR-width and MM-mean. On the other hand, (D-Ala²)-GIP significantly reduced TWCmean and CSI-width. None of the other compositional parameters were affected by (D-Ala²)-GIP administration. (D-Ala²)-GIP-Tag significantly augmented XLR-mean, XLR-width and MM-mean in the bone matrix, and reduced CSI-width. However, (D-Ala²)-GIP-Tag treatment failed to modify tissue water content. Multiple regression analyses highlighted the importance of modifications of collagen maturity (mean and heterogeneity) in the changes in ultimate load and work-to-fracture (Table 1).

DISCUSSION

Bone fragility is recognized as a comorbidity of type 2 diabetes mellitus, which has been suggested to improve with better glycemic control. However, as some glucose-lowering therapies have previously been shown to increase fracture risk [START_REF] Schwartz | Diabetes, bone and glucose-lowering agents: clinical outcomes[END_REF], it is a prerequisite to assess the safety and potential beneficial effects of such pharmacological interventions. GIP is an established incretin hormone, but it also makes an important contribution to bone strength [START_REF] Gaudin-Audrain | Glucose-dependent insulinotropic polypeptide receptor deficiency leads to modifications of trabecular bone volume and quality in mice[END_REF][START_REF] Mieczkowska | Glucose-dependent insulinotropic polypeptide (GIP) receptor deletion leads to reduced bone strength and quality[END_REF]. In the present manuscript we have evaluated the effects of GIP analogue therapies on bone strength in a rodent model of diet-induced obesity. We compared two GIP analogues: the native peptide with a D-alanine in position 2, as well as this stable GIP analogue with an additional C-terminal modification of 6 aspartic acid residues, conferring a bone tropism and called GIP-Tag [START_REF] Mabilleau | Efficacy of targeting bonespecific GIP receptor in ovariectomy-induced bone loss[END_REF]. As evident from the present results, both peptides were capable of increasing bone matrix mineralization and maturity in vitro, but also reduced osteoclast activation in the presence of high glucose concentration (equivalent to 25 mmol/L), suggesting that these peptides might present with interesting properties to prevent fragility fracture in diabetes. These observations are in agreement with previous published literature showing a reduction of bone resorption following GIP infusion in humans, which seems independent of glycemia [START_REF] Christensen | Glucose-Dependent Insulinotropic Polypeptide (GIP) Inhibits Bone Resorption Independently of Insulin and Glycemia[END_REF]. Interestingly, the addition of the aspartic residues Tag at the C-terminus of GIP did not appear to affect its biological activity. Previously, we have shown GIP-Tag to exhibit similar binding activities at the GIP receptor as (D-Ala²)-GIP or GIP1-42 [START_REF] Mabilleau | Efficacy of targeting bonespecific GIP receptor in ovariectomy-induced bone loss[END_REF]. Furthermore, activation of intracellular pathways such as adenylate cyclase, intracellular calcium, p38α, CREB, AMPKα2 and STAT2 have been evidenced with parent and GIP-Tag peptides [START_REF] Mabilleau | Efficacy of targeting bonespecific GIP receptor in ovariectomy-induced bone loss[END_REF].

In the present study, administration of (D-Ala²)-GIP and (D-Ala²)-GIP-Tag to obese pre-diabetic animals, resulted in higher ultimate bending load and work-to-fracture with no modifications of cortical (and trabecular) bone microarchitecture. This suggests that improved bone strength was due to modifications of bone tissue composition. This is very interesting when compared with the effects of exenatide, a known GLP-1 mimetic that not only improved mechanical strength but also induced changes in cortical bone microarchitecture and bone tissue composition [START_REF] Mansur | The GLP-1 Receptor Agonist Exenatide Ameliorates Bone Composition and Tissue Material Properties in High Fat Fed Diabetic Mice[END_REF].

Multiple regression analysis suggested that ultimate load and work-to-fracture were associated with changes in collagen maturity, but none of other changes in bone tissue composition. This observation shows the limit of studying tissue composition by Fourier transform infrared imaging. Thus, although it is indicative of changes in several bone composition parameters, changes in bone composition do not ultimately result in positive modifications of bone strength at the whole bone level. Blood glucose levels were decreased at the end of the study by both GIP analogues. With respect to the incretin action of GIP, it is not surprising to observe such effects with (D-Ala²)-GIP, although assessment of circulating insulin concentrations would be required to confirm this. However, the effect of (D-Ala²)-GIP-Tag that rapidly homes to bone, is perhaps surprising and suggests the possibility that some bone-produced factors are responsible for the control of glycemia. In this regard, development of specific assays to assess the pharmacokinetic profile of (D-Ala²)-GIP and (D-Ala²)-GIP-Tag would be interesting, but beyond the scope of the current study. In addition, it is worth noting that glucose-induced insulin secretion was not altered by either of the treatments. One related factor could be osteocalcin, which has been shown to function as an insulin secretagogue [START_REF] Lee | Endocrine regulation of energy metabolism by the skeleton[END_REF][START_REF] Ferron | Osteocalcin differentially regulates beta cell and adipocyte gene expression and affects the development of metabolic diseases in wild-type mice[END_REF]. However, previous reports suggested that GIP reduced osteocalcin expression and translation [START_REF] Kainuma | Attenuation by incretins of thyroid hormone-stimulated osteocalcin synthesis in osteoblasts[END_REF]. This area will require further investigation in the future. Whether a GIP-osteocalcin reciprocal relationship exists, as observed with GLP-1 [START_REF] Mizokami | Oral administration of osteocalcin improves glucose utilization by stimulating glucagon-like peptide-1 secretion[END_REF][START_REF] Mizokami | Osteocalcin induces release of glucagon-like peptide-1 and thereby stimulates insulin secretion in mice[END_REF], will also require more detailed study.

The similarity between the effects of (D-Ala²)-GIP and (D-Ala²)-GIP-Tag suggest that skeletal GIP receptors are required for beneficial effects of GIP on bone properties. This is in contrast to what was seen in ovariectomy-induced bone loss and osteopenia as (D-Ala²)-GIP-Tag was far less effective than (D-Ala²)-GIP in enhancing bone strength and quality [START_REF] Mabilleau | Efficacy of targeting bonespecific GIP receptor in ovariectomy-induced bone loss[END_REF]. However, another possible scenario would be represented by the amelioration of the glycemic profile with both peptides, and it is possible that the beneficial effects on bone are not represented by direct action of GIP on bone parameters but rather by reduction of the negative effects of chronic hyperglycemia. Further studies are required to fully ascertain the mechanism behind the beneficial bone effects. In addition, the impact of both treatments of GIP secretion, and ambient endogenous GIP levels, would be interesting to consider in follow-up investigations.

In conclusion, this study has demonstrated the potential of GIP analogues not just in the treatment of diet-induced obesity and prediabetes, but also to improve bone strength through effects on bone tissue composition. Further studies in humans are required to ascertain whether GIP analogues might represent a viable option in preventing bone fragility in diabetes.

FIGURE LEGENDS

TABLES

Figure 1 [

 1 Figure 1 [On-line color only]: Effects of (D-Ala²)-GIP or (D-Ala²)-GIP-Tag on bone cells in vitro. MC3T3-E1 cells were differentiated in the presence of ascorbic acid, Na2HPO4 and a final glucose concentration of 25 mM in the presence of vehicle, (D-Ala²)-GIP or (D-Ala²)-GIP-Tag. (A) Mineralization extent (B) lysyl oxidase activity and (C) enzymatic collagen crosslinking in the extracellular matrix. Peripheral blood mononuclear cells were also differentiated in osteoclasts in the presence of vehicle, (D-Ala²)-GIP or (D-Ala²)-GIP-Tag. (D) Number of newly-formed osteoclasts/well and (E) extent of surface resorption by differentiated osteoclasts. Vehicle: white bars, (D-Ala²)-GIP: black bars, (D-Ala²)-GIP-Tag: grey bars. *: p<0.05 vs. vehicle-treated cultures. Data are presented as median ± interquartile range from 3-5 independent experiments.

Figure 2 [

 2 Figure 2 [On-line color only]: Metabolic parameters at the beginning and end of the 42 days period of pharmacological intervention. (A) Body weight and (B) blood glucose were assessed three days before administering vehicle or peptides and then every three days during the 42-day period of the study. (C) Blood glucose was measured prior to and after i.p. administration of glucose (18 mmol/kg) at t = 0 min on day 43 in fasted animals and (D) incremental area under the glucose curves (AUC) between 0 and 105 min was computed. *: p<0.05 vs. saline-treated animals. NS: not significant. Data are presented as median ± interquartile range for n=6/group.

Figure 3 [

 3 Figure 3 [On-line color only]: Effects of (D-Ala²)-GIP or (D-Ala²)-GIP-Tag on bone mechanical resistance. Biomechanical response of femurs was tested by three-point bending. (A) Ultimate load, (B) Stiffness, (C) Post-yield displacement (PYD) and (D) Work-to-fracture. Saline: white bars, (D-Ala²)-GIP: black bars, (D-Ala²)-GIP-Tag: grey bars. *: p<0.05 vs. saline-treated animals. Data are presented as median ± interquartile range for n=6/group.

Figure 4 [

 4 Figure 4 [On-line color only]: Effects of (D-Ala²)-GIP or (D-Ala²)-GIP-Tag on bone structural properties. Tibias were scanned by high resolution X-ray microcomputed tomography and trabecular and cortical bone microarchitecture were evaluated in the proximal metaphysis. (A) Three-dimensional reconstruction of proximal metaphysis, (B) trabecular and (C) cortical bone parameters. BV/TV: trabecular bone volume density, Tb.N: trabecular number, Tb.Th: trabecular thickness, Tb.Sp: trabecular separation, Tt.Ar: total tissue area, Ma.Ar: Marrow area, Ct.Ar: cortical bone area, Ct.Th:

Figure 5 [

 5 Figure 5 [On-line color only]: Effects of (D-Ala²)-GIP and (D-Ala²)-GIP-Tag on bone matrix composition. The composition of the mineral and organic phase of the bone matrix was evaluated by Fourier transform infrared imaging. The of each parameter over the bone frequency area was computed and curve fitted with a Gaussian curve. The average and full width at half maximum of the Gaussian distribution were also calculated for each parameter and presented as "mean" and "width". XLR: collagen maturity, XST: mineral crystallinity, CSI: crystal size index, MM: mineral-to-matrix ratio, C/P: carbonate content, Ac.P: acid phosphate content, TWC: tissue water content, CG: collagen glycation. Saline: white bars, (D-Ala²)-GIP: black bars, (D-Ala²)-GIP-Tag: grey bars. *: p<0.05 vs. salinetreated animals. Data are presented as median ± interquartile range for n=6/group.

  Figure 1

Table 1 : Multiple regression analyses between bone compositional and mechanical parameters.

 1

	Dependent variable Model adjusted R² Model p value	Parameter Coefficient	p value
	Ultimate load	0.846	<0.001	Intercept	-45.6	0.015
				XLR_mean	25.1	0.022
	Work-to-fracture	0.408	0.033	XLR_width	14.435	0.033

peak, corresponding to the C=O vibration of pMMA, and the amide I band, located between 1590-1710

ACKNOWLEDGEMENTS

The authors are grateful to Nadine Gaborit and Stéphanie Lemière (University of Angers, GEROM-LHEA, Institut de Biologie en Santé, Angers, France) for their help with microCT. This work was supported by grants from the Irish Endocrine Society as well as University of Ulster Research Challenge Fund and Proof of Principle Funding Programs.

fasting blood glucose • (D-Ala²)-GIP and (D-Ala²)-GIP-Tag enhanced ultimate load and work-to-fracture • (D-Ala²)-GIP and (D-Ala²)


Fabrizio Sossan

Enrica Scolari

Rahul Gupta

Mario Paolone

Solar Irradiance Estimations for Modeling the Variability of Photovoltaic Generation and Assessing Violations of Grid Constraints: a Comparison between Satellite and Pyranometers Measurements with Load Flow Simulations

Global horizontal irradiance (GHI) is typically used to model the potential of distributed photovoltaic (PV) generation. On the one hand, satellite estimations are non-pervasive and already available from commercial providers, but they have a limited spatiotemporal resolution. On the other hand, local estimations, e.g., from pyranometers, sky-cameras and monitored PV plants, capture local irradiance patterns and dynamics, but they require insitu monitoring infrastructure and upgrading the asset of electrical operators. Considering that in most power systems, PV generation is typically the aggregated contribution of many distributed plants, are local GHI estimations necessary to characterize the variability of the power flow at the grid connection point (GCP) and detect violations of the limits of voltages and line currents accurately? To reply, we consider GHI measurements from a dense network of pyranometers (used to model the ground truth GHI potential), satellite estimations for the same area, and information about a medium and low voltage distribution system. We perform load flows at different levels of installed PV capacity and compare the nodal voltages, line currents, and the power at the GCP when the irradiance is from pyranometers and when from satellite estimations, deriving conclusions on the necessity, or not, of highly spatiotemporally resolved irradiance estimations.

I. INTRODUCTION

Current operational practices of electrical distribution utilities related to PV generation typically involve the use of solar irradiance estimations to evaluate its production potential. The kind of estimations usually depends on the level of aggregation. Distribution system operators (DSOs) typically use model-based estimates of PV generation for feasibility studies when connecting new PV power plants to verify that the grid can withstand the power injections without determining violations of statutory voltage limits and cables' ampacities (see e.g. [START_REF] Ricciardi | Defining customer export limits in pv-rich low voltage networks[END_REF] ). Irradiance forecasts might be relevant for DSOs to schedule the operation of utility-scale storage facilities to mitigate the impact of PV generation on their grids, or if specific grid codes (see, e.g., the notion of balance group in the Swiss grid code 2 ) to improve the forecast of their aggregated demand when they are penalized for incurring in unbalances.

At the level of the transmission system operators (TSOs) and load balance responsible, the need for forecasts of stochastic generation and demand stems from the requirements of dispatching conventional power plants and allocating power reserves in day-ahead and intra-day markets, typically at 1-hour or sub-hour resolution; due to the low time resolution of the dispatch processes, irradiance estimations from numerical weather predictions a) Electronic mail: fabrizio.sossan@mines-paristech.fr.; Formerly at DESL, EPFL.

(NWPs) fed by satellite information may be used (e.g., 3 ).

However, future increased generation capacity levels from renewables and reduced values of system inertia will determine more significant reserve requirements and the compression of the time scale at which ancillary services need to be delivered. Future grid codes will promote the development of markets for enhanced regulation and congestion management (e.g. [START_REF]AEMO Australian Energy market Operator[END_REF][START_REF] Hong | Fast frequency response for effective frequency control in power systems with low inertia[END_REF] ), fostering the participation of distributed energy resources, including renewable generation, in the provision of grid ancillary services. Since PV generation is the fastest-growing among stochastic renewables in distribution grids, it is relevant to wonder whether satellite estimations will be able to capture PV production patterns accurately or if local estimations should complement them.

To reply to this question, we consider a low voltage (LV) and a medium voltage (MV) distribution feeder, and GHI estimations from both a dense network of pyranometers and satellites overlaid on a common geographical meshgrid. With load flow simulations, we determine nodal voltages, line currents, and power flows in the grids, first, when GHI is from pyranometers measurements (i.e., ground-truth case) and, second, with satellite estimations and compare the results. The comparison is performed from two standpoints: an analysis on the nodal voltages and lines currents to verify whether they respect the prescribed limits (i.e., a critical operational concern for distribution system operators), and evaluation of the volatility levels of the power flow at the grids' connection points (i.e., relevant in the context of scheduling power reserves).

The use of local and satellite estimations of GHI was considered in the existing literature in application to power system studies, however little effort was devoted to comparing their performance and to understand whether higher or lower spatiotemporal resolution plays a role. For instance, Ruf et al. [START_REF] Ruf | Quantifying residential pv feed-in power in low voltage grids based on satellite-derived irradiance data with application to power flow calculations[END_REF] compare the total power production of a cluster of PV plants spread on a 0.6×0.7 km area against HelioClim-3 v.5 (HC3v5 ) estimates and ground observations from a single weather station at 11 km distance; they conclude that satellite estimations lead to a better performance than distant local observations. Leloux et al. [START_REF] Leloux | Review of the performance of residential pv systems in france[END_REF] used satellite GHI estimates to characterize the performance of PV plants, with no focus however on grid voltage and current constraints and without considering local irradiance measurements. Conversely, Nguyen et al. [START_REF] Nguyen | High pv penetration impacts on five local distribution networks using high resolution solar resource assessment with sky imager and quasi-steady state distribution system simulations[END_REF] discusses applying local GHI measurements to determine the PV hosting capacity of MV grids and the position of on-load tap changer transformers for voltage regulation, with no performance comparison with satellite estimations. Rikos et al. [START_REF] Rikos | Stability and power quality issues in microgrids under weather disturbances[END_REF] investigates the use of satellite estimations to characterize the variability of PV generation and anticipate shortcoming of generation, concluding that preemptive control based on early warnings from satellites are possible if below critical levels of installed PV capacity and ramping rates of PV generation, with no performance comparison against local GHI estimations. Lave et al. [START_REF] Lave | High-frequency irradiance fluctuations and geographic smoothing[END_REF] propose a model to characterize the geographical smoothing of distributed irradiance time series based on measurements from 6 pyranometers within 3 km distance. They observed that the variance of the spatially averaged series is six times smaller than single series for time scales up to 256 sec; above that, the spatial smoothing effect decreases due to series becoming more correlated.

The main contribution of this paper is a quantitative comparison of the performance of satellite versus local estimations of GHI for detecting grid-constraint violations and estimating the variability at the GCP. As specific numerical findings depend on the topology and physical characteristics of the lines, it is generally not possible to derive analytic principles to determine apriori whether the use of readily available satellite estimations is reliable. However, the proposed results provide useful indications to determine whether local estimations should be considered.

The remainder of this paper is organized as follows. Section II describes the modelling principles and GHI data sources, Section III details the metric adopted for the comparison, sections IV and V presents and discusses the results, and Section VI outlines the conclusions.

II. SIMULATION MODELS AND DATA

A. A common geographical reference for measurements and networks

To the best of our knowledge, an experimental setup with geolocated information on both the electrical grid infrastructure and local irradiance is not publicly available at the time of this study. Thus, we collect information from different sources and overlay them on a com-mon geographical mesh grid, as shown in Fig. 1. This allows to establish the link between GHI and grid topologies (as explained in Section II E), and to calculate the PV generation at each node. There are 3 layers of information: electrical networks, GHI pyranometer measurements, and GHI satellite estimations. The geographical area of interest extends for approximately 48 km 2 . A mesh grid with 0.040 km resolution (this approximates the smallest distance in the model, i.e. the 0.035 km long segments of the LV grid) is overlaid on the area and used to compute a uniform map of approximated GHI with nearest-neighbor (NN) interpolation. The time resolution of the time series is 30 s, chosen as a trade-off between capturing fast GHI fluctuations, and tractable computational times.

The specifications of the benchmark electrical systems (described in Section II B) do not report the position of the nodes but only the lines length; these are used to establish the relative positions of the nodes with respect to each other; once they are determined, the location of the grids as a whole is chosen so that each falls in the perimeter delimited by the outermost pyranometers so to guarantee reliable GHI interpolations. FIG. 2. Locations and labels of the pyranometers, MV nodes (in the left plot), and LV nodes (right). Excluded pyranometers refer to units that occasionally underestimated irradiance and were disregarded, as explained in the text. The distance between nodes {2, 3}, {7, 8}, and {4, 5} is less than the meshgrid resolution and share the same location on the meshgrid.

tions for the European grid benchmarks for three-phase MV and LV distribution systems [START_REF]Benchmark systems for network integration of renewable and distributed energy resources[END_REF] (figures 4 and 5). The MV grid is a 14-node 20 kV/25 MVA meshed system; the LV grid is an 18-node 0.4 kV/400 kVA radial system. The power flow at the GCP, nodal voltages and line currents of each grid are computed with single-phase load flows as a function of the grid admittance matrix and nodal injections. We thereby assume balanced three-phase systems: the inclusion of unbalances in the system is not of particular interest as most of the PV installations are interfaced to the grid with three-phase power converters. The nodal injection at each node is the algebraic sum of the electrical demand and PV generation, determined as discussed in the next two paragraphs.

We assume that excess PV generation flows in the upper-level grid. PV self-consumption schemes are not considered because not relevant in the context of comparing local versus satellite estimations.

b. Demand profiles In-line with the CIGRE specifications [START_REF]Benchmark systems for network integration of renewable and distributed energy resources[END_REF] , the active power demand at the nodes is modeled by scaling the time series shown in Fig. 3 for the nominal demand of each node (reported in tables I and II) and according to the type of load (i.e., domestic or commercial). The same re-scaled series is applied to the demand at all nodes by considering perfect correlation among them. This choice is because our focus is comparing the effects of using satellite estimations against using local irradiance, and modeling spatial correlation would introduce a new element of variability in the analysis. Reactive power injections are modeled using a constant power factor, reported in Table I.

c. PV generation CIGRE specifications 11 include a moderate quantity of connected PV generation capacity. To reproduce a future scenario with large installed capacity levels of PV generation, we apply a procedure similar to work [START_REF] Sossan | Large scale deployment of pv units in existing distribution networks: Optimization of the installation layout[END_REF] to determine the maximum PV installed power than the grid can host without violating grid constraints (i.e., so-called PV hosting capacity of the grid). We consider fixed (non-tracking) PV plants. Nominal capacity, and tilt and azimuth per node are specified in tables I and II . Plants' tilt and azimuth are, for the MV grid, the same for all the plants and optimized to achieve the highest yearly capacity factor, as typical for largescale installations. For the LV grid, they are randomly distributed around the optimal configuration to account that in LV feeders PV generation is typically from nonuniform rooftop installations. The output of each PV power plant is determined as a function of its installed capacity, tilt and azimuth, and local GHI as detailed in Section II D. PV plants are assumed operated at unitary power factor.

In addition to the case where the grid is at the limit of its PV hosting capacity (that we denote as Scenario 33), we consider a more extreme second case where the PV installed capacity at each node is scaled by factor 4/3 to analyze the sensitivity of the results with respect to the amount of PV generation (Scenario 43).

C. Data sources for the global horizontal irradiance

We describe here the irradiance estimations and measurements processing. We introduce the following notation: t = 1, 2, . . . , T denotes the index of the T time inter- nometer and satellite estimations for time t and location p are respectively denoted by x tp and xtp .

Pyranometers measurements and processing

a. Data selection Pyranometers data are from the HOPE measurements campaign [START_REF] Madhavan | Shortwave surface radiation budget network for observing small-scale cloud inhomogeneity fields[END_REF] and refer to the period from 2 April to 24 July 2013. The setup included 99 openfield installations with silicon photodiode pyranometers distributed over an area of approximately 48 km 2 . Measurements were labeled with quality flags according to the cleanliness of the sensor and horizontal alignment. We consider measurements marked as good and low incident angles with solar zenith angle smaller than 70 • to avoid inaccurate readings. To filter out potential other bad data, we remove GHI values within a ±30 s range from a missing or non-valid data point. Additionally, we verify the presence of pyranometers with outlier measurements with respect to others and exclude them as explained next. Let x (d) p be the reading of pyranometer p averaged over day d, formally:

x (d) p = 1 |D d | t∈D d x tp , (1)
where |D d | denotes the number of elements in D d , and xd the average over space:

xd = 1 99 99 p=1 x (d) p . (2)
Daily average series x p , ∀p, d is shown in Fig. 6. It was observed that pyranometers number 38, 43, 58, 60, 71, 72 (shown in Fig. 2) occasionally underestimated the irradiance with respect to others; setting up an empirical symmetrical acceptance range 0.75-1.25 (shown with the dashed lines in Fig. 6 and that contains 99.5% of the measurements) was enough to exclude them. Variations around 1 in Fig. 6 can be explained by local weather patterns and non-uniform cloud passing over the area. To avoid periods with spatially sparse measurements, we consider time intervals where at least 70% of pyranometers are available. b. Data processing GHI, originally at 1 Hz, is downsampled to 30 s by decimation (samples average would smooth quick variations) to reduce the number of load flows to perform, and NN interpolation is applied to compute one data-point per pixel of the mesh grid. It is worth noting that, on the one hand, temporal down-sampling reduces temporal variations, thus biasing the analysis in favor of satellite estimations. On the other-hand, spatial NN interpolation exacerbates spatial correlation, thus biasing the results in favor of pyranometer estimations.

c. Accuracy of measurements Madhavan et al. [START_REF] Madhavan | Shortwave surface radiation budget network for observing small-scale cloud inhomogeneity fields[END_REF] characterize the accuracy of pyranometers readings and report uncertainty levels of ±4.4 and 35.5 Wm -2 at 50 and 1000 Wm -2 , respectively, and zenith angle smaller than 70 • , i.e., ±8% and 3.5%. We anticipate that, as it will be shown in Section IV C, the impact of measurements uncertainties on the performance assessment is negligible.

Satellite estimations

a. Data selection and processing We use commercial HC3v5 GHI estimations from Soda Service. They are obtained by estimating the clear-sky index using 15-min resolution satellite images from second generation Meteosat, applied to McClear's clear-sky model including information on water vapor, aerosols, and ozone levels. The commercial provider achieves a higher time resolution with an interpolation procedure [START_REF] Lefèvre | Mcclear: a new model estimating downwelling solar radiation at ground level in clear-sky conditions[END_REF] that considers the position of the sun in the sky, but not cloud movements. We query data at a 1-min resolution for each pyranome-ter location. We achieve the target 30-sec resolution with linear interpolation.

b. De-biasing Vernay et al. [START_REF] Vernay | Characterizing measurements campaigns for an innovative calibration approach of the global horizontal irradiation estimated by helioclim-3[END_REF] correct the bias with a non-linear regression model fitted on a year of groundbased GHI measurements. This procedure is not applicable here because pyranometer data are for 4 months. Alternatively, we correct the bias by subtracting, for each pyranometer location, the difference between the satellite estimation for that location and the ground-based measurements, averaged on a daily basis, as explain next. The error of satellite estimation is:

e tp = xtp -x tp (3)
The daily bias is the average of the error over one day:

e (d) p = 1 |D d | t∈D d e tp . (4)
Satellite estimations are de-biased by applying:

x pt = x pt -e (y(t)) p , (5)
where y(t) denotes the index d of the day to which time interval t belong (y(t) = d | t ∈ D d ). Fig. 7 shows the empirical probability distribution function of the error of debiased satellite estimations: variations around zero reflect inaccuracies due to, e.g., reduced spatiotemporal resolution and parallax, as for instance reported in work [START_REF] Marie-Joseph | On the applicability of the heliosat-2 method to assess surface solar irradiance in the intertropical convergence zone, french guiana[END_REF] . Fig. 8 shows the standard deviation of the error as a function of the pyranometer location: the variability is within the range 125-145 Wm -2 . It indicates that the variations of the satellite estimation error is nearly uniform over space. This can be explained by the movements of clouds that determine similar shading patterns on the area over time and comparable levels of variability across locations. The de-biasing process (3)-( 5) might be inaccurate for late morning and evening hours, when irradiance values are lower than for the rest of the day. However, irradiance values with zenith angle larger than 70 degrees are removed. Moreover, in the context of the proposed analyses on variability and violations of grid constraints, low irradiance values are not expected to play a relevant role.

It is also worth noting that the process (3)-( 5) is not viable in real-time applications because it is non-causal (it entails calculating the bias before knowing the realizations) and requires a dense network of measurements, that is generally unavailable. However, as the comparison focuses on quantifying the effect of spatial resolution and not of bias, we consider ideally de-biased satellite estimations.

c. Quality of satellite estimations Qu et al. [START_REF] Qu | Improving helioclim-3 estimates of surface solar irradiance using the mcclear clear-sky model and recent advances in atmosphere composition[END_REF] compare HC3v5 estimations against ground-based measurements from 23 sites and report relative bias and rRMSE in the ranges -4/4% and 14/38%. The same metrics (defined in work [START_REF] Qu | Improving helioclim-3 estimates of surface solar irradiance using the mcclear clear-sky model and recent advances in atmosphere composition[END_REF] ) calculated for the current case study score in the ranges -7/3% and 28/41% at 15 min resolution, in-line with results reported in the literature.

D. Computation of the nodal PV injections

We model PV generation with the empirical model [START_REF] Skoplaki | On the temperature dependence of photovoltaic module electrical performance: A review of efficiency/power correlations[END_REF][START_REF] Sossan | Unsupervised disaggregation of photo-voltaic production from aggregated power flow measurements of heterogeneous prosumers[END_REF][START_REF] Nespoli | An unsupervised method for estimating the global horizontal irradiance from photovoltaic power measurements[END_REF] that accounts for the temperature effect by scaling the value of the POA irradiance. Let b denote the index of the electrical nodes, t the time interval, I ( b) t in Wm -2 the POA irradiance. The PV generation and the cell temperature of the PV plant at node b and time t are modeled by (assuming constant unitary efficiency with respect to the irradiance 21 ):

P V b t = I ( b) 1000 1 + γ T (b) cell,t -25 • P (b) nom (6)
T (b) cell,t = T (b) air,t + βI ( b) t , (7)
where

P (b)
nom is the nominal capacity at the standard test condition (STC), T air,t is the air temperature from measurements, and β = 0.038 and γ = -0.0043 are plant-specific parameters which were obtained by averaging the values for the close roof mount and open rack configurations [START_REF] Kratochvil | Photovoltaic array performance model[END_REF] and values for polycrystalline modules 23 , respectively.

The POA irradiance

I ( b) t
is computed by applying a transposition model as a function of the GHI, plant tilt/azimuth and time of the day as described in II E 0 b. In turn, the GHI at the location b of the plant is computed as described next.

E. Computation of plane-of-array irradiance

This step allows us to establish the link between plant position and georeferenced measurements of the GHI and to determine the POA irradiance of each PV plant. For each plant, we determine first the geographical coordinates in the mesh grid that encompass its full extension as a function of the PV capacity. These coordinates are used to query GHI with the nearest neighbor interpolation scheme. Finally, GHI is transposed to get the POA irradiance. This whole process is detailed in the following.

a. Computation of GHI We first retrieve the coordinates of the pixels encompassing the full extension of the PV plant as a function of its nominal power capacity. For convenience, we assume square PV plants and centred on the pixel of the electrical node to which they are connected. The pixels which include the plants' extension are identified starting from the plant location by performing consecutive iteration of morphological dilatation until the pixels span an equal or larger area than the plant footprint area (see Fig. 9). Assuming horizontal panels, the footprint area A b (m 2 ) of a plant with capacity P (b) nom is approximated by:

A b = P (b) nom /d, (8)
where d = 0.15 is the PV conversion efficiency. The plant area in number of pixels approximated to the nearest larger integer is:

n b = ceil A b /r 2 , ( 9
)
where r is the resolution of the meshgrid in meters. Let (i b , j b ) be the coordinates of bus b, k a non-negative integer denoting the algorithm iteration, S k . The set of pixel coordinates that contains the PV plant at node b (i.e., gray pixels in Fig. 9) is found by:

S (b) 0 = {(i b , j b )} (10) Repeat for k=1,2,. . . S (b) k = S (b) k-1 ⊕ 1 3×3 (11)
Until

|S (b) k | ≥ n b . ( 12
)
Once the set is found, it is used to query the GHI values at each coordinate following the nearest neighbor interpolation scheme. GHI values are finally averaged to get equivalent value for the whole plant. Partial shading effect will be considered in future works where more dense measurements will be available. b. Computation of POA irradiance For transposing irradiance starting from the GHI, we adopt the modeling tool-chain implemented in PVlib [START_REF] Stein | Pvlib: Open source photovoltaic performance modeling functions for matlab and python[END_REF] to allow reproducibility of the results. For the computation of the direct normal irradiance from GHI, we use the DISC model [START_REF] Maxwell | A quasi-physical model for converting hourly global horizontal to direct normal insolation[END_REF] . For transposition as a function of the panels tilt and azimuth of each plant, we use the anisotropic Hay-Davies model [START_REF] Hay | Calculation of the solar radiation incident on an inclined surface[END_REF] , which is known to have better performance than isotropic models [START_REF] Hay | Calculating solar radiation for inclined surfaces: Practical approaches[END_REF] . Results were also tested for the Perez model [START_REF] Perez | Modeling daylight availability and irradiance components from direct and global irradiance[END_REF] ; the change of the transposition model (that is x -coordinate y -coordinate FIG. 9. A geographical mesh-grid with the PV plant location (green pixel), plant boundary (green line) and set of pixels encompassing the plant extension (gray and light gray pixels), obtained, in this case, with 2 iterations of morphological dilation.

known to be a critical component, see e.g. work [START_REF] Lave | Evaluation of global horizontal irradiance to plane-of-array irradiance models at locations across the united states[END_REF] ) was however found to impact results marginally, with slight changes on false negative detections and unaltered levels of variability-underestimation factors.

III. METRICS FOR THE COMPARISON

The objective of the comparative analysis is to evaluate whether the use of satellite estimations leads to underestimate reserve capacity requirements and if they can detect violations of grid constraints reliably.

For the former objective, we rely on the notion of variability of the power export, and we verify whether satellite estimations can record, or not, the same levels of variability as pyranometers. For the latter, we verify whether satellite estimations can correctly detect violations of grid constraints, and we analyze durations and inter-arrival times of false detections. The metrics are described in the following.

a. Variability-underestimation factor of satellite estimations Let x be a time series (for instance, irradiance) with elements x 1 , x 2 , . . . and x its moving average on a 15-minute interval. We define the variability of x as the square root of the variance (Var) of the difference of the two series

variability(x) = Var (x -x), ( 13
)
whereis element-wise difference among vectors. In words, we say that (13) measures the intensity of the variations of x around a trend. Let x pyr and x sat refer to estimates from pyranometers and satellite, respectively. The variability-underestimation factor is the ratio between them:

variability-underestimation factor(x) = = variability(x pyr ) variability(x sat ) . ( 14
)
To further exemplify, we could assume that the variability is the metric that a TSO/load balance responsible uses to allocate power reserves. For instance, the higher the variability is, the larger the reserves should be. Assuming that the TSO uses satellite-derived information to allocate reserves capacity, the variability underestimation factor ( 14) can be interpreted as a scale coefficient that denotes the "inadequacy" of the allocated reserves with respect to real requirements. The variability-underestimation factor is calculated for the active power flow at the grid substation transformer (denoted by Power at the GCP), i.e., the quantity that a TSO would use to allocate reserve capacity. In addition, for the sake of following analysis, it is computed for the following time series:

• Spatially averaged GHI: GHI at the PV plants' locations averaged over all locations. It allows to compare the quality of satellite estimations compared to those of pyranometers;

• Spatially averaged POA: POA irradiance on PV plants' panels averaged over all locations. It allows to verify how the transposition model affects variability;

• Total PV production: the sum over all the plants of the PV power output. It allows to verify how local irradiance patterns are transformed in PV generation accounting for the modeled spatial distribution of PV panels;

b. False positive and negative detections of grid constraints violations A voltage violation occurs when the voltage magnitude at any node of the network is outside the range 0.9-1.1 pu. As we consider single-phase equivalent load flows and PV generation (i.e., voltage magnitudes monotonically increase with active power injections), violations occur in the form of trespassing the upper-bound 1.1 pu; we refer to this as over-voltage. A violation of a cable's ampacity limit occurs when its current magnitude exceeds 1 pu; we refer to this as overcurrent. For assessing the quality of satellite estimations in detecting over-currents and over-voltages, we introduce the notion of false negative and false positive. A false negative is when satellite estimations do not lead to a constraint violation and the pyranometers estimations do. A false positive is when satellite estimations lead to a constraint violation and the pyranometers estimations do not. False negatives are normalized by dividing for the total number of true positives, and false positives for true negatives. From a grid operator's perspective, false negatives are more of a concern than false positives because over-currents eventually lead to destroying cables.

c. Duration and inter-arrivals of false negative events It is of interest to characterize the duration of groups of false negatives and the inter-arrival between them to understand whether false detections tend to be sporadic or instead persistent over specific time intervals. The notion of duration and inter-arrival is graphically introduced in Fig. 10: the duration d 1 , d 2 , . . . is the series with the number of elements in each group of false negatives, whereas the inter-arrival i 1 , i 2 , . . . is the distance between each false negative and the next.

IV. RESULTS ON A 3-DAY DATA SET

We analyze three days of measurements to identify and illustrate some mainstream results, that will be then further verified in Section V on a longer data set. Measurements were randomly selected among the available daily time series to reflect characteristic irradiance patterns because, as evidenced by the analysis, this was found to impact on results. In particular, they refer to a nearly clear-sky day (Day A, April 22), a partly cloudy day with frequent GHI variations (Day B, May 12), and a second nearly clear-sky day (Day C, July 2), with mild GHI variations in the afternoon. GHI series averaged over all the locations are shown in Fig. 11, where the thick lines and the envelopes show respectively averages and min/max values of the estimations over space. b. False-positive and false-negative detections Table IV summarizes the number of over-currents and overvoltages in the ground-truth load flows with pyranometers estimations, and of false detections with satellite estimations. Scenario 33 features no over-voltage and, in day B, a few over-currents, which are all undetected by satellite estimations, as denoted by 100% over-current false negatives; false detections in Day B are because satellite estimations do not record the large irradiance peaks that lead to over-currents. In Scenario 33, satellite estimations do not score any false-positive detection.

A. Medium voltage voltage network a. Variability-underestimation factor

Scenario 43 features higher PV installed capacity, and few over-voltages and more over-currents than Scenario 33. Days B and C are still the most critical for false negatives due to highly volatile irradiance patterns. Scenario 43 reports less over-current false negatives than Scenario 33 because the larger PV capacity determines higher average line currents; thus smaller irradiance peaks are enough to trigger a valid detection of a grid constraint violation.

Scenario 43 reports some false positives, denoting that satellite estimations overestimates the PV production potential in some instances. This is the case, e.g., in the central part of Day A. as visible in Fig. 11, due to residual positive bias. c. Irradiance and over-current patterns during false negatives We perform this analysis based on Scenario 43 and partly-cloudy Day B because it reports the largest number of over-current false negatives.

The upper-panel plot in Fig. 12 shows the CDFs of the GHI estimation error (3) of satellite estimations during correct, false negative and false positive detections. It shows that the GHI error for false negatives and false positives is respectively smaller and larger than for true detections. This in-line with the previous insight, according to which irradiance underestimations lead to false negatives and overestimations to false positives.

The bottom-panel plot in Fig. 12 shows the CDFs of the estimation errors of line currents for satellite estimations (i.e., the difference between line currents estimated from satellite and pyranometers). Current errors are in per unit of the ampacity limit of the respective cable, e.g. 0.5 pu corresponds to half of the cable rating. The plot denotes that negative errors (current underestimations) are conducive to false negatives, and positive errors to false positives, as to be expected. It also denotes that the error magnitude is generally larger for false negatives than for false positives. Nearly 20% of false-negative episodes have current errors larger than 0.6 pu. Fig. 13 shows the histograms of the duration (in the first row) and of the inter-arrival times (second row) of over-current false negatives. We first discuss the duration by referring to the histograms of the first row. The duration of events with persistent false negatives is generally less than 10 samples (i.e., 5 minutes), except for a single event in nearly clear-sky day C, which is 20-sample long. Duration patterns are similar in the 3 days; thus, they do not seem to exhibit specific variations according to the sky conditions.

The histograms in the second row of Fig. 13 show interarrival times mostly smaller than 10 samples, denoting that false positives tend to be in the proximity of each other. Day C features longer inter-arrivals, up to 120 samples, denoting that false negatives occur in 2 groups that are 1 hour apart.

B. Low voltage network

The results for the LV network were found similar to the MV case, except for some differences ascribed to the reduced spatial scale. In this section, we briefly summarize the findings and focus on details only in the case of major differences.

a. Variability-underestimation factor Table V reports the variability-underestimation factor for days A, B and C in Scenario 33. Compared to the MV case, the LV grid scores larger variability-underestimation due to its smaller extension and reduced spatial smoothing. Largest underestimation factors happen in the nearly clear-sky day. Average transmission losses (not reported in the table) are approx. 8% of the total active power flow at the GCP. TABLE V. LV grid: variability-underestimation factors. Metric Total PV production (no T) refers to PV production modeled without considering temperature effects.

Variability-underestimation

Day A Day B Day C factor of Spatially averaged GHI 12.6 3.8 13.5

Spatially averaged POA 13.5 3.9 13.4

Total PV production 9.0 3.9 14.3

Power at the GCP 7.4 3.9 13.7

Total PV production (no T) 9.0 4.1 14.8 b. False positive and false negatives Table VI reports the number of over-currents and over-voltages in the load flows with pyranometers, and of false detections with satellite estimations. Similarly to the MV case, days B and C are the most critical for false-negative detections. False positives occur only in the scenario with the largest PV installed capacity. The CDFs of the satellite estimations error and line currents error confirm identical findings as before and they are omitted.

c. Duration and inter-arrival times Statistics of duration and inter-arrival of over-current false negatives were found comparable with the previous case, denoting that differences in the geographical dimension at this scale does not seem to play a role. Similar results are also confirmed by the analysis of the longer data set in the next section.

C. Impact of pyranometers' inaccuracies on results

To assess if the accuracy of pyranometers' measurements could impact on results, we have repeated the simulations of the MV grid by altering all irradiance value by ±5%(i.e., given by averaging the accuracy figures reported in Section II C). It was verified that results were minimally impacted, and conclusions unchanged.

V. RESULTS ON 25 DAYS OF DATA

We consider 25 days of measurements and verify if considerations developed in the former section apply. Data are shown in Fig. 14, along with their classification according to prevailing sky conditions (clear-sky, nearly clear-sky, partly cloudy, cloudy), performed by visual inspection.

a. Variability Fig. 15 shows the variabilityunderestimation factor of Total PV production against the variability of spatially averaged GHI of pyranometers (defined in Section III). They denote that, i), variability is generally higher in partly cloudy and nearly clear-sky days than in cloudy and clear-sky days, and ii), largest to lowest underestimations are in nearly clear-sky, partly cloudy, clear-sky and cloudy days. Results for the MV and LV grids are similar, however the former features smaller axis values than the latter due to its larger extension and spatial smoothing. In this regard, it is interesting to observe that despite the MV grid uses information from 11 pyranometers on a surface of 15 km 2 (whereas the LV grid uses 2, which are less than 400 meters far apart, and extends for 0.035 km 2 ), the variability of its spatially averaged GHI is only 30% smaller than for the LV grid, denoting an overall small spatial smoothing effect.

The variability-underestimation factor was defined in (14) as the ratio of the variability recorded by pyranometers; hence, it expresses the proportion among the two. To complement that information, we define the variability difference variability difference(x) = (15) variability(x pyr )variability(x sat ), (16) which can be interpreted as the magnitude of the uncaptured variations. Fig. 16 shows the variability difference of Total PV production against the variability of spatially averaged GHI of pyranometers for the MV network. They stand in a nearly linear relationship, denoting that, the higher the variability of the pyranometers is, the more satellite estimations underestimates the tives due to large spatial variations of irradiance measurements, as the orange envelope in Fig. 14 shows, that are not captured by satellite estimations. Other cloudy days are not shown in the plots because either score zero false positive and negative, or their irradiance patterns do not determine grid violations. Similar considerations apply to the LV grid (right panel), with the differences that the proportion of FPs is lower and more cloudy days scores false negatives, denoting higher susceptibility to local irradiance conditions.

c. Duration and inter-arrival times The analysis on the duration and inter-arrival on the 25-day data set did not show specific patterns with respect to sky-conditions and are not reported here. Also, differences between the MV and LV grids did not emerge in the comparison.

d. Congested lines

The plot in the upper left panel of Fig. 18 shows the occurrences of over-currents in the MV grid according to sky conditions. Over-currents mostly happen in nearly clear sky days (where irradiance feature very large values due to cloud enhancement) followed by clear-sky days and partly cloudy days. Congested lines are the number 3 and 15-17, which, from Fig. 4, are the closest to the GCP and the most critical because they are subject to the aggregated power flows. The plot in the bottom left panel of Fig. 18 shows the occurrences of over-currents false negatives according to sky conditions. They occur mostly in nearly clear-sky and partly cloudy days, and less in clear-sky days, when satellite estimations are able to detect irradiance dynamics better.

The plots in the right panels of Fig. 18 show the same results for the LV network. Congested lines are 1-3 and 11-16, which are those close the GCP and in the lateral feeders that interface a large amount of generation, respectively. Distribution of both over-currents and false negatives show the same features as discussed for the MV case.

VI. CONCLUSIONS AND PERSPECTIVES

a. Summary of the work We investigated the performance of HelioClim-3 v.5 satellite estimates of GHI compared to ground-based measurements from pyranometers (HOPE campaign) in applications to modeling the variability of distributed PV generation and detection of grid constraints violations (nodal voltages and line currents). The analysis was performed with load flow simulations considering the grid and demand specifications of the CI-GRE benchmarks for LV and MV systems. Distributed PV systems were sited and sized with a procedure from the literature to saturate the hosting capacity of the grids. PV generation at each node was modeled as a function of installed capacity and GHI with a transposition model. Satellite estimations were de-biased on a daily basis to retain the focus the analyses on variability only. Pyranometers measurements were post-processed to remove outliers and potential other bad data close to missing measurements. Analysis were conducted at a 30 sec resolution.

b. Key results The analyses were carried out for 28 daily profiles of the irradiance with several typical sky conditions: clear-sky, nearly clear-sky, partly cloudy, and cloudy. It was found that:

• at the edge of the grids' PV hosting capacity, satellite estimations could not detect mild over-current episodes in both LV and MV grids, leading to false negative estimations;

• false negatives were more likely in partly cloudy (when local cloud enhancement phenomena affect PV production patterns) and cloudy days. False positives were more likely in nearly and clear-sky days due to a residual local positive bias of satellite estimations;

• false negative detections were persistent for long time intervals (e.g., tens of minutes) due to the residual negative bias of satellite estimations, as well as sporadic (e.g., single false negatives with inter-arrival times of minutes) due to fast fluctuations of the irradiance that were not captured by satellite estimations;

• satellite estimations underestimated the variability of the total PV production and the power flow at the GCP in both LV and MV grids. Underestimations were less severe in the MV grid case due to the larger geographical scale. Non-uniformly installed PV generation capacity and temperature were found affecting variability patterns, suggesting that these phenomena should be taken into considerations when, for instance, plan reserve requirements.

c. Implications and perspectives

There is no disputing on the practical value of satellite estimations. However, results suggest that, when using irradiance estimations with similar attributes as those considered (e.g., spatiotemporal resolution and parallax), proper counter measurements or data reconciliation strategies with local estimates should be adopted to limit the shortcomings described above. These considerations also extend to forecasts, that are typically elaborated starting from satellite images and are delivered on relatively coarse geographical mesh grids. Implications are case specific. For TSOs and balance responsible parties, coarse irradiance estimates in confluence with high levels of installed PV generation capacity can lead to underestimating the variability of PV generation, thus underestimating reserve requirements, which should be therefore scheduled adopting conservative margins (this can also apply to DSOs if grid codes require them to compensate deviations from pre-defined dispatch plans). DSOs, when planning the connection of new PV facilities, should also implement conservative margins from grid constraints or adopt realistic models of spatiotemporal variations, especially when near the PV hosting capacity of the grid, to avoid false negative detections. In the context of real-time and predictive control of grids for congestion management, dispatch, voltage control, and peak shaving, results suggest that spatiotemporal coarse estimations could be complemented with realtime local information at a higher resolution to perform better decisions.

Current results were elaborated from load flow simulations. The uncertainty related to pyranometers measurements and two different anisotropic transposition models (Hay-Davies and Perez) were found to impact on numerical results marginally without altering conclusions. If field data from an experimental setup with many gridconnected distributed PV plants will become available, a comparison with the proposed results would be of definite interest to confirm findings as they would inherently include all modeling uncertainties (i.e., direct/diffuse separation, transposition, panels' spectral response and efficiency, and distribution grid's parameters). Also, the effects of uncorrected bias should be further investigated; even if this is not expected to play a significant role on variability, it could certainly have a dominant impact on false positive and false negative detections of grid constraints violations.

FIG. 1 .

 1 FIG. 1. Location of pyranometers and electrical networks. The zoomed region shows the LV grid. Information are from different sources and overlaid on a common geographical mesh grid.

  FIG.3. Active power profiles of the demand in per unit (pu) of the nominal value of the demand.

FIG. 4 .FIG. 5 .

 45 FIG.4. The CIGRE benchmark grids for MV systems[START_REF]Benchmark systems for network integration of renewable and distributed energy resources[END_REF] ; breaker S1 is operated closed.

p

  /x d , that essentially tell us how a measurements of a pyranometer compare with respect to average. The cumulative distribution function (CDF) of x(d)

FIG. 6 .

 6 FIG.6. CDF of normalized daily irradiance recorded by the pyranometers after filtering outliers.

FIG. 7 .

 7 FIG. 7. Histogram of satellite estimation errors.

2 )FIG. 8 .

 28 FIG.8. Standard deviation (SD) of the satellite estimations error over time as a function of the location. Missing points refer to pyranometers that were excluded.

  k the set of coordinates at k, ⊕ the Minkowski sum, 1 3×3 the all-one 3×3 matrix , and |S (b) k | the cardinality, or number of elements, of S (b)

1 FIG. 10 .

 110 FIG.10. Duration (d1, d2, d3) and inter-arrival (i1, i2, i3) of false detections.

1

  FIG. 12. MV grid: CDFs of GHI errors of satellite estimations (upper panel) and of line current errors of satellite estimations (lower panel) in Day B and scenario 43.

  FIG.13. Duration (first row) and inter-arrival (second row) of false negative events for the MV grid case.

)FIG. 14 .

 14 FIG. 14. Fifteen non-consecutive daily scenarios of GHI: D, E, F, G, H (20, 23, 24, 25, and 26 April 2013), I, J, K, L, M, N, P (2, 4, 5, 6, 12, 15, and 16 May 2013), Q, and R (4 June and 3 July). Values with solar zenith angle larger than 70 • are omitted. Solid lines denote spatial averages and color-filled areas denote the envelope between maximum and minimum values over space.

FIG. 15 .FIG. 16 .

 1516 FIG.15. Variability underestimation factor of total PV generation of satellite estimations vs variability of Spatially averaged GHI of pyranometers for the MV (left) and LV grid (right). Alphabetic labels refer to the 25 days of data.

FIG. 17 .

 17 FIG. 17. False-negative versus false-positive detections of over-currents of satellite estimations for the MV (left) and LV grid (right). Alphabetic labels refer to the 25 days of data.

FIG. 18 .

 18 FIG. 18. Top panel: congested LV lines according to sky conditions in the MV (left) and LV grid (right). Bottom panel: LV with false negative detections according to sky conditions in the MV (left) and LV grid (right).

TABLE I .

 I MV grid: nominal loads, power factor (pf ), PV installed capacity at STC, and ID and distance of the nearest pyranometer per node.

	Node		Demand		PV Generation		Nearest
		Residential Commercial	Power	pyranometer
		MW	pf	MW	pf	MWp	ID dst (m)
	1	15.3 0.98 5.1	0.95	5.71	10	700
	2	0	-	0		5.02	49	769
	3	0.28 0.97 0.26 0.85	4.39	77	256
	4	0.44 0.97	0	-	2.23	86	126
	5	0.75 0.97	0	-	1.28	26	165
	6	0.56 0.97	0	-	0.09	86	268
	7	0	-	0.09 0.85	0	78	544
	8	0.6 0.97	0	-	0	78	641
	9	0	-	0.67 0.85	1.13	89	400
	10	0.49 0.97 0.08 0.85	1.62	91	200
	11	0.34 0.97	0	-	5.71	91	89
	12	15.3 0.98 5.28 0.95	5.71	55	408
	13	0	-	0.04 0.85	0	59	322
	14	0.22 0.97 0.39 0.85	5.02	18	544
	Total 34.28		11.91		37.91

TABLE II

 II

	. LV grid: nominal load, power factor (pf ), installed
	capacity at STC and configuration (i.e., azimuth and tilt) of
	PV generation, and ID and distance of the nearest pyranome-
	ter per node.
	Node	Demand	PV Generation			Nearest
		Residential Power Azimuth Tilt pyranometer
		kW	pf	kWp	•	•	ID dst (m)
	1	200 0.95	227	167	34 75	80
	11	15	0.95	114	161	43 75	165
	15	52	0.95	107	185	38 75	165
	16	55	0.95	112	182	42 86	113
	17	35	0.95	26	190	41 86	80
	Total 404		586

Table III

 III power at the GCP reflects the same trends as former quantities, however with different values than Total PV production because, i ), the time-varying demand profile impacts on the variability and, ii ), of path-dependant transmission losses. In Scenario 33, grid losses are approximately 5% of the total power flow in average.

	Scenario 43 denoted increased PV power production lev-
	els. Results are not reported because they do not bring
	new relevant insights on variability.
	re-
	ports the variability-underestimation factor for days A,
	B and C in scenario 33. It can be observed that:
	1. largest underestimations happen in nearly clear-sky
	days A and C, denoting that satellite estimations
	cannot spot highly variable irradiance components
	caused by, e.g., small clouds, because of reduced
	spatiotemporal resolution and parallax effects;
	2. underestimation levels of spatially averaged
	POA irradiance are different than of GHI. This is
	to be expected because irradiance transposition is
	nonlinear;
	3. underestimation levels of Total PV production
	are different than of Spatially averaged POA due
	to temperature and non-uniformly distributed PV
	generation capacity that amplifies differently local

irradiance patterns; to understand the role played by these two factors, Total PV production (no T) reports the underestimation factor when the temperature effect on PV generation is not modeled. In this case, the PV production is directly proportional to the local POA irradiance. It can be seen that modeling temperature leads to reducing the underestimation factor; 4.

TABLE III .

 III MV grid: variability-underestimation factors. Metric Total PV production (no T) refers to PV production modeled without considering temperature effects.

	Variability-underestimation	Day A Day B Day C
	factor of
	Spatially averaged GHI	9.8	2.3	8.3
	Spatially averaged POA	10.9	2.2	8.3
	Total PV production	12.5	2.6	9.5
	Power at the GCP	5.6	2.6	8.2
	Total PV production (no T)	13.2	2.7	10.0

TABLE IV .

 IV MV grid: occurrences (%) of over-voltages and over-currents with pyranometers, and false positive (FPs) and false negative (FNs) detections of satellite estimations. Jul 2013). Values with solar zenith angle larger than 70 • are omitted. Solid lines denote spatial averages and color-filled areas denote the envelope between maximum and minimum values over space.

			Pyranometers max/min	Satellite max/min	Satellite spatial average	Pyranometers spatial average
	)
	GHI (W/m 2	0 500 1 000	A 8h	A 12h	A 16h	B 8h	B 12h	B 16h	C 8h			C 12h	C 16h
						UT Time (day label and hour, UTC)
	FIG. 11. Satellite (SAT) and pyranometer (PYR) estimations of GHI averaged over all locations for 3 days: A (22 Apr 2013),
	B (12 May 2013), and C (2
							PV penetration scenario	Scenario 33	Scenario 43
							Metric	Day		A	B	C	A	B	C
							Over-voltages (PYR)		0	0	0	0	0.32	0
							Over-currents (PYR)		0	0.03	0	6	1.57 3.54
							Over-voltage FNs (SAT)	-	-	-	-	100	-
							Over-current FNs (SAT)	-	100	-	16.6 72.4 58.4
							Over-voltage FPs (SAT)	0	0	0	0	0	0
							Over-current FPs (SAT)	0	0	0	0.4	0.2	1.7

TABLE VI .

 VI LV grid: occurrences (%) of over-voltages and over-currents with pyranometers, and false positive (FPs) and false negative (FNs) detections of satellite estimations.

	PV penetration scenario	Scenario 33		Scenario 43
	Metric	Day	A	B	C	A	B	C
	Over-voltages (PYR)	0	0	0	0	0.24 0.03
	Over-currents (PYR)	0	0.01	0	1.59 1.44 1.05
	Over-voltage FNs (SAT)	-	-	-	-	100	100
	Over-current FNs (SAT)	-	100	-	10.1 89.9 53.1
	Over-voltage FPs (SAT)	0	0	0	0	0	0
	Over-current FPs (SAT)	0	0	0	0	0	0.2

ACKNOWLEDGMENTS

The research was partly funded under the FURIES and JASM programs of the Swiss Competence Center on Energy Research. Authors warmly thank Dr. Deneke (Leibniz Institute for Tropospheric Research, Germany) for providing data of the HOPE campaign, Dr. Schmidt (EHF Department, University of Oldenburg, Germany) for having pointed out the HOPE campaign, and Prof. Jan Kleissl (University of California San Diego, USA) and the anonymous reviewers for providing us with relevant comments.