Thierry Feuillet

Etienne Edouard Cossart

Hélène Charreire

Arnaud Banos

Hugo Pilkington

Virginie Chasles

Serge Hercberg

Mathilde Touvier

Jean-Michel M Oppert

Hybridising geographically weighted regression and multilevel models: a new approach to capture contextual effects in geographical analyses

come   L'archive ouverte pluridisciplinaire

Introduction

Geospatial data analysts often face two major challenges when conducting quantitative geographical analyses. These challenges relate to the two main geographical effects that characterize all geographical data: the spatial effect (mainly controlled by distance) and the contextual (or 'platial' when contexts are geographical) effect [START_REF] Arcaya | Area variations in health: A spatial multilevel modeling approach[END_REF]Wolf et al., 2021). From a theory-driven perspective, the spatial effect relates to spatial dependence, according to which the strength of spatial interaction between two things depends on their mutual distance [START_REF] Tobler | A computer movie simulating urban growth in the Detroit region[END_REF]. The contextual or platial effect, on the other hand, is related to spatial heterogeneity -proposed by [START_REF] Goodchild | GIScience, Geography, Form, and Process[END_REF] to be the second law of Geography -since this effect is directly due to the fact that spatial contexts differ from one another, and therefore that belonging to a given spatial context implies specific consequences on objects and processes operating within it. From a statistical viewpoint, spatial dependence is usually linked to spatial autocorrelation, while spatial heterogeneity leads to spatial nonstationarity, i.e., to means, variances and covariances that vary from a place to another.

The demarcation line between spatial effects and contextual effects is not always a clear one, but a fundamental distinction between the two effects can be drawn depending on how we consider geographical space: either as continuous (thus revealing distance effects), or as discrete (thus capturing contextual effects). This dichotomy echoes the distinction between space and place discussed by human geographers for many years [START_REF] Kearns | Space in its place: Developing the link in medical geography[END_REF][START_REF] Roche | Geographic information science II: Less space, more places in smart cities[END_REF], Wolf et al., 2021). The place perspective uses location to form groups, while the spatial perspective looks at proximity between observations. As raised by [START_REF] Arcaya | Area variations in health: A spatial multilevel modeling approach[END_REF], the respective merits of these two approaches are rarely compared for one given dataset.

This dichotomy is also to be found in statistical techniques that are commonly used to model geographical data. Two major families of statistical modelling frameworks are used by most studies dealing with geographical effects. The first encompasses, within a place-based perspective, all the mixed models that include spatial context considered as a random effect (the so-called hierarchical or multilevel models, [START_REF] Goldstein | Multilevel statistical models[END_REF]. The second comprises, within a space-based perspective, spatially explicit models including spatially varying coefficient models, i.e., those considering geographical coordinates for each observation to estimate relationships, such as geographically weighted regressions (GWR, [START_REF] Brunsdon | Geographically weighted regression: A method for exploring spatial nonstationarity[END_REF] upon which we will focus here.

From a theoretical viewpoint, each of these families of models has strengths and limitations.

Multilevel models, popular in the education, social and health sciences since the 1990s [START_REF] Chaix | Comparison of a spatial perspective with the multilevel analytical approach in neighborhood studies: The case of mental and behavioral disorders due to psychoactive substance use in Malmö, Sweden[END_REF][START_REF] Duncan | Context, composition and heterogeneity: Using multilevel models in health research[END_REF], provide through their hierarchical structure a sound inferential framework for modelling contextual effects while simultaneously controlling for individual effects. Moreover, random slope models allow the capture of heterogeneity in direction or intensity of relationships between places [START_REF] Jones | People, places and regions: Exploring the use of multi-level modelling in the analysis of electoral data[END_REF]. Thus, these multilevel models are able to avoid both ecological and atomistic fallacies, as outlined by [START_REF] Owen | Under examination: Multilevel models, geography and health research[END_REF]. A major pitfall of these models is that the boundaries of spatial contexts must be established before the estimation of the model, and these boundaries determine in turn potential contextual effects. In the social sciences, administrative boundaries are commonly used to define spatial contexts and thus to group observations. Contextual effects can thus be over-or underestimated if administrative boundaries do not match with the actual and relevant geographical contexts related to the outcome under study [START_REF] Merlo | A brief conceptual tutorial of multilevel analysis in social epidemiology: Linking the statistical concept of clustering to the idea of contextual phenomenon[END_REF]. According to [START_REF] Owen | Under examination: Multilevel models, geography and health research[END_REF], this problem refers to spatial design. Many scholars, particularly in the field of neighbourhood effect studies, have outlined this spatial mismatch between administrative neighbourhood boundaries used by default as spatial contexts because of ease of access, and the true causally relevant context [START_REF] Diez-Roux | Neighborhoods and health[END_REF][START_REF] Petrović | Freedom from the tyranny of neighbourhood: Rethinking sociospatial context effects[END_REF]. According to the modifiable areal unit problem (MAUP), such a spatial mismatch may have considerable consequences on subsequent statistical associations, including in multivariate models [START_REF] Fotheringham | The modifiable areal unit problem in multivariate statistical analysis[END_REF], leading to seek out the "optimal" zoning prior to statistical modelling [START_REF] Openshaw | An empirical study of some zone-design criteria[END_REF].

On the other hand, the increasingly popular GWR constitutes a powerful solution to explore spatial nonstationarity without any preset geographical contexts. Based on local regressions, GWR provides maps of regression coefficients that can be further used to identify spatial contexts with specific and sometimes opposing trends in relationships [START_REF] Fotheringham | Geographically weighted regression: The analysis of spatially varying relationships[END_REF]. However, this method does not allow us to account for the aggregation of observations within places through variance component decomposition, like multilevel models do. Therefore, it ignores the effect of "belonging together" among those observations that fall within geographic boundaries (Figure 1).

To overcome the limitations of these methods, while retaining their strengths, this study investigates the implementation of a methodological workflow combining the two approaches into a two-step procedure. The key idea is to (i) delineate GWR-based spatial contexts, (ii) incorporate them into a multilevel modelling framework as a random effect, and (iii) compare the resulting model, in terms of statistical quality and fit, with the one based on administrative boundaries. These GWR-based spatial contexts are hypothesized to capture "true" causal spatial contexts (i.e., those that actually do affect human behaviour but that remain unknown) better than administrative boundaries, within a data-driven perspective (Fotheringham & Sachdeva, 2022). Some studies have previously explored varying degrees of hybridization between spatial and multilevel models, but not in the to the extent we describe here. For instance, [START_REF] Chen | Using multilevel modeling and geographically weighted regression to identify spatial variations in the relationship between place-level disadvantages and obesity in Taiwan[END_REF] studied socio-ecological determinants of obesity in Taiwan by first estimating a multilevel model to adjust mean odds ratios (ORs) by administrative units, then using these ORs as a dependent variable in a GWR. Other authors used spatial models as a first step to better define residual spatial dependence structures in multilevel models [START_REF] Dong | Spatial random slope multilevel modeling using multivariate conditional autoregressive models: A case study of subjective travel satisfaction in Beijing[END_REF][START_REF] Janko | Extending multilevel spatial models to include spatially varying coefficients[END_REF][START_REF] Park | A spatially filtered multilevel model to account for spatial dependency: Application to self-rated health status in South Korea[END_REF] or to account for spatial dependence in adding spatial autoregressive terms [START_REF] Dong | Spatial Autoregressive Models for geographically hierarchical data structures[END_REF]. Recently, [START_REF] Hu | Uncovering spatial heterogeneity in real estate prices via combined hierarchical linear model and geographically weighted regression[END_REF] proposed a combination of multilevel and GWR models (abbreviated as HLM-GWR) for modelling spatial data with identical geographical coordinates. In their study on real estate prices in China, a multilevel model was used to consider observations within the same locations as grouped, while simultaneously capturing spatially varying relationships of the group-level variables through GWR. Closer to our study, [START_REF] Arcaya | Area variations in health: A spatial multilevel modeling approach[END_REF] incorporated spatially explicit information into their multilevel models to define spatial contexts in a study on life expectancy in the US. Rather than considering counties (level-1) within states (level-2), they replaced the administrative level of US states with spatial patches solely based on proximity around counties. However, they still based their neighborhood structure scheme (level-1) on administrative units (counties).

The main originality of our article resides in the willingness to overcome pre-set geographical boundaries to construct spatial contexts. These contexts are built in a two-stage data-driven perspective, (i) First, letting GWR exhibit spatial non-stationarity, using it to draw contexts independently of administrative boundaries and (ii) including the resulting outputs in a multilevel model as a grouping variable.

This methodological procedure was applied to the geographical variation in a wellestablished indicator of obesity using data from participants of the French web-based Nutrinet-Santé cohort dataset residing in the Paris region. In this cohort, each participant was located at the residential address. Obesity is recognized as an important public health issue leading to increased morbidity and mortality (WHO Europe Obesity Report, 2022) and the study of obesity distribution fits well with a contextual analysis including both space and place effects. Indeed, a large body of research has highlighted that the prevalence of obesity reveals spatial patterns that differ according to scale [START_REF] Swinburn | Dissecting obesogenic environments: The development and application of a framework for identifying and prioritizing environmental interventions for obesity[END_REF], and that both direct (e.g. individual obesity-related behaviors such as diet and physical activity) and indirect (e.g. built and social environment) obesity drivers are spatially non-stationary [START_REF] Chen | Using multilevel modeling and geographically weighted regression to identify spatial variations in the relationship between place-level disadvantages and obesity in Taiwan[END_REF][START_REF] Feuillet | Spatial heterogeneity of the relationships between environmental characteristics and active commuting: Towards a locally varying social ecological model[END_REF][START_REF] Feuillet | Influence of the urban context on the relationship between neighbourhood deprivation and obesity[END_REF][START_REF] Oshan | Targeting the spatial context of obesity determinants via multiscale geographically weighted regression[END_REF].

In the following sections, we first briefly describe the overall methodological workflow, before describing each step in more detail, as well as the data we use, and finally reporting and discussing the main results and possible extensions of the proposed methodology.

Overall methodological workflow

The methodological workflow is based on three main steps (Figure 2): (i) First a GWR model is estimated; (ii) A spatially constrained multivariate clustering is then applied to the GWR coefficient map in order to delineate data-driven spatial contexts; Finally (iii) a multilevel model including these contexts as a random effect is estimated, and then compared to a counterpart model incorporating administrative units as a grouping variable.

Dataset and variables

In this study we use data from the ongoing French Nutrinet-Santé web-based cohort [START_REF] Hercberg | The Nutrinet-Santé Study: A web-based prospective study on the relationship between nutrition and health and determinants of dietary patterns and nutritional status[END_REF], restricting to participants residing in the Paris region. Launched in 2009, this cohort aims to provide information about relationships between nutrition and health among more than 100,000 participants aged 18 years or older, who completed through a secured website a set of questionnaires assessing their socioeconomic and health-related characteristics. Residential addresses were obtained from all participants, geocoded to the parcel or street levels. We focused on the area covering Paris inner city and its three surrounding "départements" called the "Petite Couronne" (Figure 3), 6.9 million inhabitants (2022 French census), i.e., a population density of ~8000 hab./km² (> 20,000 hab./km² in Paris inner city). Data for 9,086 participants were available for analysis in the study area after the removal of missing residential addresses. The outcome variable is the body mass index (BMI), calculated as weight [kilograms]/height² [meters]. BMI is the indicator most commonly used in epidemiological and surveillance settings to assess body weight and excess adiposity defining obesity (WHO Europe Obesity Report, 2022). Although BMI is simply a function of weight and height it is strongly linked to ill health. Given that this article focuses on a methodological dimension, a minimal set of explanatory variables have been included, in order to emphasize computation and complexity (we will address the question of less parsimonious models in the discussion). Therefore, in the first step of the GWR, only two individual variables were considered: sex (categorical) and age (continuous and scaled beforehand) of participants. While socioeconomic profiles are known to be strongly associated with obesity, we did not consider them since they are typically spatially patterned.

The linear functional form between age and BMI was visually checked using a spline-based generalized additive model.

Step 1 -Geographically weighted regression

We hypothesize that some predictors of a dependent variable have greater impact in some places than in others. GWR addresses this underlying assumption by fitting local regression models to each individual location [START_REF] Brunsdon | Geographically weighted regression: A method for exploring spatial nonstationarity[END_REF][START_REF] Fotheringham | Geographically weighted regression: The analysis of spatially varying relationships[END_REF]. Each local model uses an inverse distance weighting scheme such as bisquare or Gaussian kernel functions, assuming that closer observations from the calibration point are more influential than those further away. A basic GWR model extends traditional regression as follows:

(Equation 1) Where is the dependent variable at location i, denotes the geographical coordinates of the i th location (i.e. individual), is the k th independent variable at location i, is the intercept at location i, is the local regression coefficient for the k th independent variable at location i, and is a normally distributed error term at location i. The estimation of the local regression coefficients is given by the following equation (in a matrix form):

(Equation 2)

Where is the sampling matrix of independent variables and is the diagonal spatial weight matrix defining the weights given to each neighbour of an observation during the local calibration:

(Equation 3)

Where n is the number of neighbours at the i th location. Whether an observation is defined as a neighbour is given by two possible criteria: either a distance threshold, or a given number of k nearest neighbours (knn). The knn approach is advised when observations are irregular over space, to avoid too few observations in some local kernels [START_REF] Feuillet | Spatial heterogeneity of the relationships between environmental characteristics and active commuting: Towards a locally varying social ecological model[END_REF].

Subsequently, this method was used in this study. Thus, different kernel functions can be used to define neighbouring spatial weight matrices. We tested the five following kernel functions [START_REF] Gollini | GWmodel: An R package for exploring spatial heterogeneity using geographically weighted models[END_REF]:

The bi-square function was found to be the best one to minimize the GWR AIC and was therefore kept in further analyses (note that the methods and results of the GWR-related parameter optimization -knn, but also the number of clusters -will be addressed in a dedicated final method section, because these are related to the overall workflow). Finally, each set of parameters was tested for spatial non-stationnarity using the F3 test described in [START_REF] Leung | Statistical tests for spatial nonstationarity based on the geographically weighted regression model[END_REF] and implemented in the GWmodel R package [START_REF] Gollini | GWmodel: An R package for exploring spatial heterogeneity using geographically weighted models[END_REF][START_REF] Lu | The GWmodel R package: Further topics for exploring spatial heterogeneity using geographically weighted models[END_REF]. The test statistics reflect the sample variance of the estimated values of β ik and can be approximated by a F-distribution. H0 is that all β ik are equal.

Step 2 -Spatially constrained multivariate clustering of GWR estimates

The next step consists in computing a spatially constrained multivariate clustering based on the vector of local GWR estimates. In addition to the coefficients associated to the two predictors (sex and age), local intercepts were also considered in the spatial clustering.

These intercepts reflect the average BMI when sex and age are fixed and thus give the spatial distribution of the adjusted BMI. This age-and sex-adjusted BMI is itself potentially affected by contextual effects and therefore captures many unobserved spatially structured predictors, such as socioeconomic deprivation and the density of the built environment.

Several methods exist to build spatially constrained multivariate clustering. For example, in their study of housing market regionalization based on GWR coefficients, Helbich et al.

(2013) used the SKATER (Spatial Klustering Analysis by Tree Edge Removal) algorithm, developed by [START_REF] Assunção | Efficient regionalization techniques for socio-economic geographical units using minimum spanning trees[END_REF]. However, this method is resource-intensive and would make the whole parameter optimization burdensome. To overcome this issue, we used a more efficient method proposed by [START_REF] Chavent | ClustGeo: An R package for hierarchical clustering with spatial constraints[END_REF]. This method consists in building a hierarchical ascendant classification based on two matrices. The first is a dissimilarity matrix based on the vector of the three scaled GWR coefficients while the second is a neighbour matrix (contiguity) representing the spatial constraint. Once the spatial clustering was computed, the last step was to move from the punctual classified observations to polygons.

For this purpose, we polygonised points through a Dirichlet tessellation, also known as Voronoi diagram. This regionalization procedure subdivides space in n cells (called Voronoi polygons), n being the number of points (i.e. home addresses), so that every location in a given cell is closer to its generating point than to any other. We then attributed cluster membership to each Voronoi polygon, and finally combined polygons according to clusters to form spatial contexts.

Step 3 -Multilevel (or hierarchical) modelling

The GWR-based regionalization resulting from the previous step was used to group observations in a subsequent multilevel model. Therefore we have observations i = 1,...,n clustered in groups j = 1,...,J, j representing each GWR-based spatial context. In multilevel models, we can allow groups to vary through the use of either intercepts (intercept-varying models), slopes (slope-varying models), or both (complete models). Here we estimated the following complete model, including two predictors at level-1 but no group-level predictors, expressed as follows:

(Equation 4)

Where is the response variable for an individual i in a GWR-based spatial context j, are the predictors (sex and age), is the mean intercept for all the contexts, and are the mean slopes, is a normally distributed term, and are the independent and normally distributed deviations between each context and the mean relationships.

Parameter optimization and model validation

Two essential parameters must be set during the model calibrations. The first is the GWR bandwidth (i.e., the local kernels based on knn) and the second is the number of clusters.

GWR bandwidths are commonly selected using cross-validation based on prediction accuracy. However, the purpose of this study is not to maximize the GWR prediction performance, but rather to maximize the final multilevel model derived from GWR estimate clustering. Consequently, we selected optimal parameters (bandwidth and number of clusters) in minimizing the Akaike information criteria (AIC) of the final multilevel model.

Since the AIC function has potentially several local minima, it requires the implementation of a constrained optimization of a non-linear, multivariate and complex objective function. We used the generalized simulated annealing algorithm [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF], as being able to tackle these constraints and to converge more quickly than other similar methods [START_REF] Xiang | Generalized simulated annealing[END_REF]. Beyond this optimization, robustness analyses were performed on the number of clusters (from 20 to 60) and type of administrative units (EPCI and municipalities) and will be presented in the results section. Note that the administrative-based model was based on a French specific intercommunality administrative structure called Etablissement public de coopération intercommunale (EPCI, k = 31). We also compared results using municipalities (known as communes) as a grouping variable (j = 143). All the analyses have been written in R language (R Core Team, 2023). The code core structure is shared on the following GitHub page: https://github.com/tfeuillet/gean2023, but note that given our analyses involved personal health data, the code cannot be be directly reproducible.

Results

The mean age of the participants was 42.4 years and 76.8% were women. This overrepresentation of women may reflect a volunteer bias, knowing that it is established that women are more likely to participate in research studies (Andreeva 2016, Galea, 2007) and are also more health-and nutrition-conscious than men (Barebring, 2020). . Mean BMI was 23.4 kg/m2 (descriptive statistics are presented in Table 1).

GWR model outputs

GWR scaled coefficient surfaces are presented in Figure 4 and their descriptive statistics in Table 2. Coefficients are interpolated (inverse distance weighting method) for visualization purposes, and the spatial distribution of participants is mapped using kernel densities (Fig. 3A). The optimization algorithm yielded an optimal GWR bandwidth equaling 67 nearest neighbours of participants. The mean distance between each participant and its 67 th nearest neighbour is about 1 km. While the OLS model reveals that BMI is on average positively and significantly correlated with age and with being a female, GWR shows significant spatial nonstationarity for the intercepts and age coefficients, with regular opposite signs according to location (Figure 4 B,C,D). There is a strong difference between the median GWR coefficient for age, and the OLS coefficient (0.98 vs. 0.07, respectively, see Table 2), indicating that many local coefficients are higher than the mean OLS coefficient. The intercept map also displays contrasted values throughout the study site. Regarding the GWR-based regionalization, the results provide 60 spatial contexts, mapped in Figure 5, alongside the 31 administrative units used in the comparison model.

Multilevel model estimation

Results of the two multilevel models (administrative-based and GWR-based, respectively, the municipality level being discussed later) are provided in Table 3 and Figure 5. The first finding is that the AIC of the administrative-based model is higher than the GWR-based model (52326 vs. 52156), indicating a better quality of the hybrid model. The second finding is that the intraclass correlation coefficient (ICC) is more than twice as high in the GWRbased model as that for the other (5.9% vs. 2.8%), reflecting higher correlations among observations in the GWR spatial contexts than in the administrative units. Finally, the last interesting observation is the difference in slope variances (noted τ11) between the two models, equaling 0.032 in the first and 0.212 in the GWR-based model. This difference is illustrated in Figure 5 for the partial effect of age. This means that heterogeneity of relationships is stronger in GWR contexts than in administrative ones. Table 3. Results of the multilevel models. Left: model using administrative units (N = 31) as a grouping variable (i.e. random effect). Right: model using GWR-based spatial contexts as a grouping variable (N = 60). τ00 refers to the intercept variance, and τ11 to the slope variance (for age). Graphs display the relationships between age and BMI in the two models. The violin plots (bottom) show that both intercept and slope variances are higher in the GWR-based model.

Administrative-based multilevel model GWR-based multilevel model

Model robustness under alternative parameters

In multilevel modelling, random effect variances depend on the number of units in the grouping variable. To check for such a possible bias, we (i) conducted robustness analyses by varying the number of GWR clusters from 20 to 60 and (ii) estimated a multilevel model based on municipalities (j = 143) instead of intercommunalities (EPCI) as used previously.

Results showed that AIC of GWR-based models remained lower than administrative-based models notwithstanding the number of clusters (Figure 6). Note that the AIC plateaued from 45 clusters. Likewise, the model based on municipalities exhibited a higher AIC than the one based on intercommunality units, as well as lower ICC and slope variance (data not shown).

Discussion

In this study we focused on developing a novel modelling approach combining two major and complementary modelling frameworks commonly used in contextual studies (GWR and multilevel models). We applied this approach to the geography of human corpulence as assessed by the BMI. When integrating GWR-based spatial contexts as contextual settings rather than administrative units (commonly used in multilevel models), we overcame the well-known spatial design issue in multilevel modelling [START_REF] Owen | Under examination: Multilevel models, geography and health research[END_REF]. We demonstrated in our sample that our hybrid approach outperforms administrative-based multilevel models as measured through decreasing AIC and account better for contextual effects through ICC and slope variance increases.

The spatial design issue in multilevel modelling refers to the mismatch between the geographical boundaries of contextual setting used in the model and the true causally relevant spatial contexts affecting the studied outcome, that remain, in fact, unknown (Diez-administrative boundaries are used by default and this choice affects subsequent inferential results about contextual effect occurrence and strength. Therefore, delineating relevant spatial contexts in multilevel models is a crucial challenge. Conceptual thinking about sociospatial context boundaries has led to interesting debates in human geography [START_REF] Feuillet | Built environment in local relation with walking: Why here and not there[END_REF][START_REF] Kwan | The uncertain geographic context problem[END_REF][START_REF] Kwan | The limits of the neighborhood effect: Contextual uncertainties in geographic, environmental health, and social science research[END_REF][START_REF] Petrović | Freedom from the tyranny of neighbourhood: Rethinking sociospatial context effects[END_REF]. Recommendations are often made to consider idiosyncratic contexts, since the way people respond to environmental influences (i.e., context effects) are highly personal [START_REF] Kwan | The limits of the neighborhood effect: Contextual uncertainties in geographic, environmental health, and social science research[END_REF]. Likewise, activity spaces are individual, and some contextual studies have used self-defined spatial contexts (Charreire et al., 2016;[START_REF] Perchoux | Conceptualization and measurement of environmental exposure in epidemiology: Accounting for activity space related to daily mobility[END_REF]. However, such an individual view of contextual effects is not without strong methodological issues, both in terms of data acquisition in large samples, and in terms of compatibility with the multilevel modelling framework.

Therefore, the solution we proposed to define spatial contexts, based on a spatial nonstationarity and a data-driven strategy, may be considered as a workable compromise between administrative units (partially disconnected from actual contextual effects), and idiosyncratic contexts (impossible to integrate into multilevel models as a random effect).

GWR-based spatial contexts seem relevant in that they can capture unobserved spatially structured factors that affect the outcome under study (here, measures of obesity such as BMI). Such multidimensional and complex factors that shape spatial contexts are typically difficult to measure quantitatively at a supra-local scale (i.e., at a scale ensuring sufficient heterogeneity), as they relate to complex and spatially heterogeneous interactions between places and people [START_REF] Feuillet | Built environment in local relation with walking: Why here and not there[END_REF][START_REF] Fotheringham | Scale, context, and heterogeneity: A spatial analytical perspective on the 2016 US presidential election[END_REF]Fotheringham & Sachdeva, 2022). According to Fotheringham & Sachdeva (2022), spatial context "is a shorthand term for the impact of the largely unmeasurable effects of location on one's actions [...] and hence a multifaceted concept incorporating the influence of local media, family, friends, and local organizations as well as notions of traditions, persistent adverse or beneficial conditions, customs, lifestyles and psychological profiles common to an area that affect social norms, which in turn affect individual behavior" (Fotheringham & Sachdeva, 2022, p. 3). In ecological studies of obesity, the complexity of obesogenic contexts leading to a high BMI (so-called obesogenic environments) including personal and socioeconomic factors, and the built environment, into a comprehensive social ecological system has been emphasized [START_REF] Swinburn | Dissecting obesogenic environments: The development and application of a framework for identifying and prioritizing environmental interventions for obesity[END_REF] and can be transposed to other fields of human-related behaviour. We argue here that local modelling techniques such as GWR can be used to efficiently reveal spatial contexts otherwise unmeasurable, and consequently be included as a grouping variable in subsequent multilevel models to relevantly capture contextual effect in a sound inferential framework. That said, it is worth noting that spatial and contextual effects are sometimes embedded, and that "unmeasurable effects of location" can actually manifest as both spatial heterogeneity and spatial dependence. This complexity leads to be cautious with the interpretations, in particular with salient effects that could require further qualitative work to reach a sound conclusion.

Caution is also crucial because of the other possible causes of GWR-derived spatial nonstationarity, leading to potential inferential bias. The first other cause raised by Fotheringham & Sachdeva (2022) is noise due to sampling variation during the subset local calibrations, inherent to kernel-based local techniques such as GWR. Such a sampling variation may result in the spatial variability of parameters, even if relationships are constant over space.

The second cause is due to the assumption of an incorrect functional form in the relationships under study. Some studies have demonstrated that nonlinear relationships modelled as linear in spatially varying coefficient models may exhibit spatial non-stationarity [START_REF] Sachdeva | Are we modelling spatially varying processes or non-linear relationships?[END_REF]. Note that a nonlinear relationship can also reflect some spatial nonstationarity. This is typically the case when a predictor is strongly spatially patterned, e.g., in following a centre-periphery gradient, such as density or socioeconomic deprivation in European cities [START_REF] Feuillet | Modelling context-specific relationships between neighbourhood socioeconomic disadvantage and private car use[END_REF]. In this study we checked the linear form of the relation between age and BMI through using generalized additive models, a highly recommended approach [START_REF] Hastie | Generalized additive models[END_REF]. However, the functional form of a relationship can be linear on average, but non-linear in specific local kernels. Taken together, these points (noise and model misspecification) constitute a pitfall that modelers must be aware of.

Another point of discussion arises from the choice to limit the number of predictors in the GWR model, as we did in this study. We deliberately selected a limited number of explanatory variables which were only individual level variables We justify this choice by the fact that including contextual variables (e.g., characteristics of the built environment expected to influence obesity-related behaviours) in GWR would require defining individual exposure to such characteristics (typically buffers around home addresses, or administrative boundaries). This would bring us back to the initial problem raised in this study, namely, how to delimit the individual exposure to a spatial context. Thus, we recommend adding contextual variables at the step of the multilevel model, i.e., after having built the GWR spatial contexts. Moreover, this solution allows us to be parsimonious in GWR and to make the regionalization procedure easier, while also avoiding the redundancy of spatial information in GWR. Indeed, including spatially structured variables in GWR makes interpretation somewhat complicated since it layers over the spatially explicit information derived from local kernels.

We tested this hybrid approach by applying it to geographical variation in BMI, the most common indicator of obesity, but it would be worth considering generalizing it to other outcomes and other fields of geographical analyses. Neighbourhood effects are well suited to such a method, since they embrace complex spatial contexts that are not easily measurable by common means. Besides multilevel models, GWR-based regionalization has already been used for housing market segmentation and further hedonic modelling [START_REF] Helbich | Data-driven regionalization of housing markets[END_REF], and we may suppose that such an approach could also be useful for some other zoning-based analyses, such as spatial interaction models. Also, potential extensions or complements could be tested, for instance in using the multiscale GWR proposed by [START_REF] Fotheringham | Multiscale Geographically Weighted Regression (MGWR)[END_REF], which is able to capture the different spatial scales at which predictors operate. The way the regionalization is done is also a parameter that could be compared and improved. For example, the spatially clustered regression recently proposed by [START_REF] Sugasawa | Spatially clustered regression[END_REF] could be an interesting way to delineate spatial contexts as well. These issues represent open avenues as areas for research, to improve and consolidate this new hybrid approach of contextual effect modelling.

Conclusion

The issue of the delimitation of a pertinent spatial context is crucial in multilevel modelling, in order to appropriately capture contextual effects. In the social sciences, administrative boundaries are often used by default and because of a relative ease of access, even if they do not match with the actual geography of the spatial processes at work. In this study, focusing on geographical distribution of the BMI as an indicator of obesity, we suggest regionalizing GWR coefficients in order to delineate unmeasurable spatial contexts assumed to affect BMI and the relationships between BMI and individual variables (age and gender).

Finally, we use these contexts as a random effect in a subsequent random slope multilevel model. We show that the GWR-based multilevel model outperformed its administrativebased counterpart, in terms of quality and contextual effect modelling. From a theoretical perspective, this hybrid procedure also provides a means to reconcile space and place, in accounting for both the spatial effect through GWR, and the place effect through multilevel modelling. When characterizing geospatial data, these two kinds of effects are rarely considered simultaneously in statistical models despite existing theoretical soundness. It can be hoped that the methodology presented here will pave the way to further extensions and improvements, in particular in parameter calibration and robustness to other data samples and outcomes, at various scales.

Figure 1 .

 1 Figure 1. Ability of common statistical models used in geographical analysis to consider spatial non-stationarity and contextual effects. The GWR-multilevel hybrid procedure proposed in this study can simultaneously handle spatial nonstationarity and contextual effects. Note that "contextual models" refers to ordinary least square (OLS) models including contextual aggregated variables as predictors.

Figure 2 .

 2 Figure 2. Overall methodological workflow followed in this study: bw is the bandwidth used in the GWR weighting scheme, k is the number of clusters, and AIC is the Akaike information criteria. The three steps marked in the figure are those that are specifically explained in the method section below.

Figure 3 .

 3 Figure 3. Location map of the study area. Black lines delineate "départements" (92 : Hautde-Seine; 93 : Seine-Saint-Denis; 94 : Val-de-Marne) while white liness define contours of intercommunality administrative structure called Etablissement public de coopération intercommunale (EPCI), used as reference units in the subsequent multilevel models.

  Validation of the final GWR-based multilevel model was achieved by (i) comparing its AIC with the one of a usual administrative-based multilevel model and (ii) comparing the intraclass correlation coefficient (ICC, which is an indicator of the correlation between two observations taken randomly within a given spatial context or administrative unit) and slope variances (indicating the ability of the model to capture contextual effects) of the two multilevel models. The model exhibiting the highest values of ICC and slope variances was considered as the one including the most relevant spatial contexts among the contexts tested in this study.

Figure 4 .

 4 Figure 4. A: Spatial distribution of the participants (N = 9089) using kernel density. B, C and D: Maps of the scaled GWR coefficients after IDW interpolation (for intercepts, gender and age, respectively). Marginal plots (in grey) summarize the mean raster values (kernel density) using zonal statistics.

Figure 5 .

 5 Figure 5. Maps of spatial contexts (EPCI and GWR-based) used in the multilevel models.

Figure 6 .

 6 Figure 6. Robustness analyses on the number of GWR-based spatial contexts. The hybrid model outperformed the EPCI-based model whatever the number of contexts in terms of AIC minimization.

Table 1 .

 1 Characteristics of the sample (N = 9089)

	N	Mean (SD) or %	Range
	Age	9089 42.4 (14.6)	[18-86]
	Gender (female)	9089 76.8 %	n/a
	BMI [kg/m 2 ]	9089 23.4 (4.5)	[12.2-76.2]

Table 2 .

 2 GWR raw coefficient summary. OLS coefficients have been added for comparison purposes.

		GWR model (non-scaled coefficients)			OLS model
							Leung's F3-test for
							spatial non-
							stationnarity (p-
		Min	1st Qu	Median 3rd Qu Max	value)	Beta	p-value
	Intercept	20.21	22.33	23.00	23.77	28.85 <0.001	20.26	<0.001
	Gender (ref =
	male)	-4.88	0.23	1.23	2.24	7.63 >0.05	1.19	<0.001
	Age	-2.09	0.55	0.98	1.44	4.77 <0.001	0.07	<0.001
	Number of
	observation
	s	9089						9089
	R 2	0.23						0.07

[START_REF] Diez-Roux | Neighborhoods and health[END_REF][START_REF] Owen | Under examination: Multilevel models, geography and health research[END_REF]. In many studies, particularly neighbourhood studies,

Acknowledgments

The authors warmly thank Pierre Barbillon and the two anonymous reviewers for their suggestions on the first version of the manuscript.


domain: phys spi

text:
Jianping Yin

F De Gregorio

K.-S Rossignol

L Rottmann

G Ceglia

Gabriel Reboul

G Barakos

G Qiao

M Muth

M Kessler

A Visingardi

M Barbarino

F Petrosino

A Zanotti

N Oberti

L Galimberti

G Bernardini

C Poggi

L Abergo

F Caccia

A Guardone

C Testa

S Zaghi

Acoustic and aerodynamic evaluation of DLR small-scale rotor configuration within GARTEUR AG 26

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Multirotor systems are increasingly applied to Vertical Take-Off and Landing (VTOL) configurations. Multirotor systems are popular for small-scale drones used for surveillance, monitoring, search, aerial surveying and photography, amongst other applications. In the context of a growing interest in the development of urban air mobility (UAM) solutions and providing answers to fundamental questions regarding the aerodynamic and acoustic characteristics of these new vehicles, the GARTEUR Action Group HC/AG-26 [START_REF] Yin | Noise Radiation and Propagation for Multirotor System Configurations[END_REF] have been established. There are 14 partners involved in the group activities coordinated by DLR. The objectives of this AG are [START_REF] Yin | Noise Radiation and Propagation for Multirotor System Configurations[END_REF] to gain knowledge in noise generation and noise propagation of multirotor systems under the influence of the installation effects and (2) to develop and validate numerical prediction methods for the noise prediction for multirotor systems. For the purpose of the code validations, the results obtained by each partner were shared among the action group. Three data bases of test data were available from DLR [START_REF] Rossignol | Investigation of Small-Scale Rotor Aeroacoustic in DLR's Acoustic Wind Tunnel Braunschweig[END_REF][3], CIRA/CUSANO [START_REF] Nargi | Evolution of flow structures in twin-rotors wakes in drones by time-Presented at 49th European Rotorcraft Forum[END_REF] [START_REF] Ceglia | Time-resolved PIV measurements of the rotor wake in a side-byside configuration[END_REF] and Polimi [START_REF] Zanotti | Aerodynamic interaction between tandem overlapping propellers in eVTOL airplane mode flight condition[END_REF]. In addition, the first common test campaigns involving PIV measurement were carried out by CIRA/DLR [START_REF] Gregorio Rossignol | Multi-rotor wake interaction characterization[END_REF], using the DLR small rotor configurations [START_REF] Rossignol | Investigation of Small-Scale Rotor Aeroacoustic in DLR's Acoustic Wind Tunnel Braunschweig[END_REF] [START_REF] Yin | Numerical Investigations on Small-scale Rotor Configurations with Validation using Acoustic Wind Tunnel Data[END_REF]. The common experiment aims at using the dedicated capabilities of the partner's wind tunnels to improve the validation database for the simulations and at the same time to validate the experimental accuracy by repeating the same experiments in several wind tunnels. This paper will focus on the numerical comparison activities with the wind tunnel results involving only DLR small rotor configurations, including isolated and rotor/rotor configurations. For the simulation of the various test cases, the numerical approaches from each partner are applied. The aerodynamic simulations necessary for the aeroacoustic predictions are conducted with a range of fidelity numerical methods, varying from lifting line to CFD. In the case of the acoustic installation effects, a scattering code can be applied or directly extract from CFD. The introduction of the experimental approach used in the acoustic/PIV test, including the rotor model, the instrumentation and the data reduction will also be presented. The methodologies applied in the numerical simulations by the partners will then be described and analyzed to point out their strengths and weaknesses. The aerodynamic and acoustic predictions will be analyzed and compared either code to code or with available test results for the different configurations considered. The acoustic results will be presented in terms of sound pressure level (SPL) directivities, spectra and sound pressure time histories.

DESCRIPTION ON THE ACOUSTIC AND PIV

TEST SETUP

Acoustic test setup and measurement in DLR AWB

The AWB [START_REF] Pott-Pollenske | Enhanced Capabilities of the Aeroacoustic Wind Tunnel Braunschweig[END_REF] is DLR's small-scale high-quality anechoic testing facility.

It is an open-jet Göttingentype wind tunnel capable of running at speeds of up to 65 m/s and optimized for noise measurements at frequencies above 250 Hz. The nozzle is 1.2 m high by 0.8 m in width. A special rig was designed to extend the capabilities of the facility to meet the requirements of simultaneous measurements of multiple rotors under static and flight conditions, e.g. Figure 1 and Figure 2.

The main objective of the selected mechanical design is to enable the investigation of the effect of flow and shaft angle on the acoustic radiation of a broad range of propeller configurations; isolated, coaxial, tandem with vertical and lateral offset, e.g. Figure 1. The dimensions of the AWB test section, allow the investigations of rotors with a diameter of up to approximately 0.4 m. The rig is designed to allow shaft angle variations in the range 𝛼 ± 30° and testing at free-stream velocities up to 𝑈 ∞ ≈ 30 m/s. These setting values are evidently valid for the current set of propellers considered herein. The whole structure of the rig is based on standard X-95 rails and carriers. This choice allows for easy changes in configuration. The rig is fixed to a rotating axle at its centre point, i.e. on the left-hand side of Figure 2 to allow variations in rotor shaft angle. The two cameras were installed side-by-side and horizontally oriented as shown in Figure 3, to increase the sizes of the measured region.

PIV test setup and measurement in CIRA

The measurement region was illuminated by a double-head Nd-Yag laser with a pulse energy of 150 mJ at a wavelength of 532 nm and a repetition rate of 10 Hz.

Figure 3 PIV set-up for tandem case in CIRA

The flow measurements provide a vector pitch of ∆𝑥 = 1.41𝑚𝑚. The error of the PIV cross-correlation procedure is evaluated as 0.1 px as first estimation (Raffel et al. [START_REF] Raffel | Particle Image Velocimetry -A Practical Guide[END_REF]). Using the current values for the optical resolution (11.3 px/mm) and the laser doublepulse delay (25-35 μs), this provides a velocity error of V of  0.25 m/s to 0.35 m/s. For a detailed description of the PIV measurement system, the reader is referred to the work of De Gregorio et al [START_REF] Gregorio Rossignol | Multi-rotor wake interaction characterization[END_REF].

For the tandem configuration, the measured region covers between -2.13<x/R<0.29 along the radial direction and the area from slightly above the propeller disk until down to little more than two radii distance z/R=-2.1 (Figure 3, Figure 5-a). For the coaxial configuration, a different camera setup was adopted considering the presence of the fairing below the propeller blade. The PIV cameras were vertically mounted side-by-side (Figure 4) assuring a measurement region ranging along the x-axis between 0.16 < x/R < 1.

Rotor and Pylon

The rotor is a commercially available one ((Xoar model PJN). For this paper, the test results from a two-blade 13x7 rotor including a rotor support (pylon), as shown in Figure 6 are chosen for validating the numerical simulations. The rotor 13x7 represents a rotor with 13 inches or 0.33m in diameter and 7 inches in pitch. The rotor, pylon planform and the form of the section profiles are demonstrated in Figure 6. The twist and chord distribution are given in The rotors are mounted on Leopard LC5065 motors coupled to YGE 205HVT speed controllers and SM300-Series 3300 W DC power supplies. This combination allows RPM up to 13000 to be reached.

For each rotor, performance data, in terms of thrust and torque, is acquired through miniature sixcomponents load cells, Modell K6D40 from ME-Meßsysteme GmbH, mounted directly underneath the propellers. Each rotor RPM is acquired through a 1/rev signal generated by a Hall-effect sensor mounted to the rotor's shaft. This signal also serves as a trigger signal for data post-processing. Two multi-propeller configurations will be considered. An overview of these is given in Figure 5 a) and b). The so-named tandem configuration is defined by a relative position between the two rotor axes of Δx=1.18 of the rotor diameter D, Δy=0 and a vertical distance of Δz=0.25D. The coaxial configuration has Δx and Δy null, while the axial distance between the propeller disks is equal to Δz=0.25 D.

DESCRIPTION OF METHODOLOGIES APPLIED IN THE NUMERICAL SIMULATIONS BY THE PARTNERS

The numerical investigations are performed by each partner applying either in-house-developed or commercial computational tools. The aerodynamic simulations necessary for the aeroacoustic predictions are conducted with a range of fidelity numerical methods, varying from lifting line to CFD.

CIRA

The CIRA aerodynamic simulations were carried out by using the medium-fidelity code RAMSYS [START_REF] Visingardi | Validation of CIRA's rotorcraft aerodynamic modelling system with DNW experimental data[END_REF], which is an unsteady, inviscid and incompressible free-wake vortex lattice boundary element methodology (BEM) solver for multi-rotor, multi-body configurations developed at CIRA. It is based on Morino's boundary integral formulation [START_REF] Morino | A General Theory of Unsteady Compressible Potential Aerodynamics[END_REF] for the solution of Laplace's equation for the velocity potential φ. The surface pressure distributions are evaluated by applying the unsteady version of Bernoulli equation, which is then integrated to provide the forces and moments on the configuration and the surrounding obstacles. A computational acceleration is obtained by applying the module for symmetrical flows and geometries implemented in the solver and the parallel execution via the OpenMP API. The ACO-FWH solver is used for computing the acoustic free-field generated by the rotor blades. It is based on the FW-H formulation [START_REF] Williams | Sound generation by turbulence and surfaces in arbitrary motion[END_REF] described in [START_REF] Casalino | An advanced time approach for acoustic analogy predictions[END_REF] [14] and [START_REF] Barbarino | A high-fidelity aeroacoustic simulation of a VTOL aircraft in an urban air mobility scenario[END_REF]. The advanced-time formulation of Farassat 1A is employed, and the linear terms (the so-called thickness and loading noise contributions) are computed through integrals both on the moving blades surface (impermeable/rigid surface formulation). The computational acceleration is obtained by a parallel execution via the MPI API. The simulation of the aeroacoustic free-field was carried out by using the aerodynamic database evaluated by RAMSYS, and consisting of the rotor blade pressure distributions.

DLR

The free-wake panel method UPM [START_REF] Ahmed | Unsteady Panel Method Calculation of Pressure Distribution on BO105 Model Rotor Blades[END_REF][17] is based on a velocity-based, indirect potential formulation using a combination of source and vortex distribution on the solid surfaces and vortex panels in the wake. Compressibility effect of the flow is considered by applying the Prandtl-Glauert correction. The blade vortex interaction (BVI) is captured thanks to the free wake model used in UPM. The validation effort is supported by CFD TAU steady simulations on selected hover test cases. The unstructured CFD code TAU is based on the solution of the Reynolds averaged Navier-Stokes equations on hybrid unstructured meshes. The solver relies on a cell vertex scheme to discretize the mass, momentum and energy fluxes [START_REF] Schwarmborn | The DLR TAU-Code: Recent Applications in Research and Industry[END_REF]. Depending on the configurations, all interactions among propellers, rotor, fuselage, wings and tail are considered. By postprocessing the unsteady aerodynamic pressure on the blade surface, the Ffowcs Williams-Hawkings [START_REF] Williams | Sound Generation by Turbulence and Surfaces in Arbitrary Motion[END_REF] (FW-H)-equation based code APSIM [START_REF] Yin | Rotor Aerodynamic and Noise under Influence of Elastic Blade Motion and Different Fuselage Modeling[END_REF] is used to compute the noise radiation and to predict the acoustic levels on a given hemisphere surface.

ONERA

The aerodynamic simulations performed by ONERA are realized with the PUMA code [START_REF] Boisard | Numerical Analysis of Rotor / Propeller aerodynamic interactions on a highspeed compound helicopter[END_REF]. PUMA (potential unsteady methods for aerodynamics) is an unsteady lifting line / free-wake solver developed at ONERA since 2013. It is built on a coupling between an aerodynamic module and a kinematic module. The aerodynamic module relies on a lifting line method with a free-wake model using the Mudry theory [START_REF] Mudry | La théorie des nappes tourbillonnaires et ses applications à l'aérodynamique instationnaire[END_REF], which describes the unsteady evolution of a wake modeled by a potential discontinuity surface. The lifting line method relies on twodimensional airfoils characteristics through lookup tables computed preliminary by CFD with the ONERA elsA1 1 code [START_REF] Cambier | The ONERA elsA CFD software: input from research and feedback from industry[END_REF]. Some blade sweep correction and dynamic stall models are added. Concerning the kinematic module, it is based on a rigid multi-body system approach using a tree-like structure with links and articulations. It enables any arbitrary motion between the different elements. In order to speed up the computation the code has been parallelized using OpenMP and the Multilevel Fast Multipole Method has been implemented.

Concerning the numerical parameters used for the computations, they are based on previous experience. The lifting line is divided in 30 radial stations using a square root distribution. A time step of 5° was used over 25 rotor revolutions over which the last 6 are used for post-processing. The computations do not account for the rotor hub or other test rig components.

The unsteady spanwise distribution of loads obtained with PUMA are used as input for the KIM code [START_REF] Prieur | Comparison of the Ffowcs Williams-Hawkings and Kirchhoff rotor noise calculations[END_REF] [START_REF] Rahier | An efficient Kirchhoff integration method for rotor noise prediction starting indifferently from subsonically or supersonically rotating meshes[END_REF] to determine the noise emission of the rotor thanks to a Ffowcs-Williams and Hawkings formulation implemented in a non-compact advanced time approach. Since only sectional forces are available and in order not to consider noise sources compact in the chord direction, the surface pressure is reconstructed over the entire blade thanks to interpolation based on the pressure distributions computed and stored during the airfoil polar computations. 1 elsA V4.2.01 Airbus-Safran-ONERA property

Polimi

The single propeller case was simulated with both a mid-fidelity and a high-fidelity aerodynamics solver, respectively DUST and SU2. DUST is an opensource software developed by POLIMI to simulate the interactional aerodynamics of unconventional rotorcraft configurations. The code, released as free software under the open-source MIT license, relies on an integral boundary element formulation of the aerodynamic problem and on a vortex particle model of the wakes [START_REF] Tugnoli | Mid-fidelity approach to aerodynamic simulations[END_REF]. SU2 is an open-source toolkit distributed by the SU2 Foundation [START_REF] Economon | SU2 An open-source suite for multiphysicsì simulation and design[END_REF], freely available and licensed under the GNU Lesser General Public License. It uses the finite volume approach to solve partial differential equations (PDE) on unstructured meshes. It solves the Unsteady Reynold-averaged Navier-Stoke (URANS) equations to analyze typical aeronautical problems that involve turbulent flows in the compressible regime. Aerodynamic results obtained with both solvers are not trimmed. The aeroacoustic signature is computed by solving Ffowcs Williams-Hawkings [START_REF] Williams | Sound Generation by Turbulence and Surfaces in Arbitrary Motion[END_REF] (FWH) equations. The surface pressure field on the propeller computed with the two solvers is provided as input to the same acoustic module [START_REF] Galimberti | Propeller Noise Prediction Capabilities within SU2[END_REF]. The pylon for both DUST/SU2 simulations is not modelled.

Roma Tre University/CNR-INM (ROMA3)

The RM3 aerodynamic and aeroacoustic analyses rely on tools developed by the Roma Tre University unit in the last twenty years and widely validated in the past in helicopter and tiltrotor configurations [START_REF] Gennaretti | Rotorcraft comprehensive code assessment for blade-vortex interaction conditions[END_REF] [START_REF] Gennaretti | Prediction of tiltrotor vibratory loads with inclusion of wing-proprotor aerodynamic interaction[END_REF]. The aerodynamic module is based on the boundary integral formulation for the velocity potential presented in [START_REF] Gennaretti | Novel boundary integral formulation for blade-vortex interaction aerodynamics of helicopter rotors[END_REF], suited for helicopter configurations where blade-vortex interaction (BVI) occurs. This formulation is fully 3D, can be applied to bodies with arbitrary shape and motion, and allows the calculation of both wake distortion and blade pressure field. It assumes the potential field to be divided into an incident field, generated by doublets over the wake portion not in contact with the trailing edge (far wake), and a scattered field, generated by sources and doublets over the body and doublets over the wake portion very close to the trailing edge (near wake). This procedure allows one to overcome the instabilities arising when the wake comes too close to or impinges on the body. Recalling the equivalence between the surface distribution of doublets and vortices, the contribution of the wake portion experiencing BVI (far wake) is expressed in terms of thick vortex (i.e., Rankine vortices) distributions. The shape of the wake can be either assigned (prescribed-wake analysis) or obtained as a part of the solution (free-wake analysis) by a time-marching integration scheme in which the wake is moved accordingly to the velocity field computed from the potential solution. Once the potential field is known, the Bernoulli theorem yields the pressure distribution on the body that, in turn, is used both to determine the aerodynamic loads and as an input to the aeroacoustic solver to predict the radiated noise. The aeroacoustic analysis is performed by a prediction tool based on the Ffowcs Williams and Hawkings equation (FWH) [START_REF] Williams | Sound Generation by Turbulence and Surfaces in Arbitrary Motion[END_REF]. The solution of the FWH equation is achieved through the boundary integral representation known as the Farassat Formulation 1A [START_REF] Farassat | Derivation of Formulations 1 and 1A of Farassat[END_REF].

Uni Stuttgart (IAG)

For high-fidelity simulation a framework consisting of FLOWer and ACCO was used at the Institute of Aerodynamics and Gas Dynamics (IAG) at the University of Stuttgart. CFD results are obtained with the block-structured code FLOWer, originally developed by DLR [START_REF] Kroll | MEGAFLOW-a numerical flow simulation system[END_REF] and continuously extended at the IAG for rotorcraft and multirotor applications [START_REF] Keßler | Progress in IAG's Rotorcraft Simulation Framework[END_REF]. Acoustic coupling was provided by IAG's FW-H solver ACCO [START_REF] Keßler | [END_REF] which uses the transient flow data provided by FLOWer as an input.

FLOWer solves the three-dimensional, compressible RANS equations and enables unsteady flow solutions (URANS). The discretization of time and space is applied separately by the method of lines. For temporal discretization a second-order dual time stepping is used [START_REF] Jameson | Time Dependent Calculations Using Multigrid, with Applications to Unsteady Flows Past Airfoils and Wings[END_REF] with a time step of 0.5° to resolve acoustic waves, while for spatial discretization a 6th high-order weighted essentially non-oscillatory (WENO) scheme by Borges [START_REF] Borges | An Improved Weighted Essentially Non-Oscillatory Scheme for Hyperbolic Conservation Laws[END_REF] is used. Furthermore, the k-omega turbulence model was applied to close the URANS equations. Using the Chimera technique separate meshes were created for the background and rotor, utilizing hanging grid nodes to reduce the numerical expense. The spatial resolution of a single rotor was achieved by 6.6 mio cells. The rotor surface was meshed with 144 cells in radial and 80 cells in chordwise direction with a refinement towards the blade tip, leading edge and trailing edge with cell sizes of less than one per mille of the radius. The rotor mesh is extruded in wall normal direction with 52 cells. The first surface cells satisfy y+ < 1 and an extrusion up to cell sizes corresponding 10% of the chord length is applied. The spatial discretization in the background mesh was based on the resolution of the first harmonic wave length, with 15 cells discretizing the wave length of the blade passing frequency (BPF). The acoustic code ACCO is an in-house code of the IAG, which uses an acoustic integration surface for the generation of sound emissions. For the integration either the physical surfaces or a permeable surface surrounding the object of interest can be used. For the performed simulations, the physical blade data of 4 full rotor revolutions are used and the integration is achieved through the physical surface of the rotor blade, which includes all tonal sound sources.

UoG (University of Glasgow, Glasgow)

Helicopter Multi-Block (HMB3) code is employed in this study. The solver can accurately predict the aerodynamic performance, acoustics of propagation, and has been widely used in the investigation of rotorcraft flows [START_REF] Antoniadis | Assessment of CFD methods against experimental flow measurements for helicopter flows[END_REF], helicopter rotor aeroelasticity [START_REF] Steininger | Numerical Analysis of HVAB and STAR Rotor Blades Using HMB3[END_REF], and missile trajectory prediction [START_REF] Babu | Aeroelastic simulations of stores in weapon bays using Detached-Eddy Simulation[END_REF]. Moreover, a good agreement when compared to experimental results in aerodynamics, acoustics and aeroelasticity of propellers, which can be seen in previous study [START_REF] Barakos | Acoustic comparison of propellers[END_REF]. Most recently, its ability to capture the interactions of multi-rotor flows and ducted propeller flows was documented [START_REF] Qiao | Parametric study of aerodynamic performance of equivalent ducted/un-ducted propellers based on high-fidelity computational fluid dynamics[END_REF]. HMB3 solves the Unsteady Reynolds Averaged Navier-Stokes(URANS) equation in integral form using the Arbitrary Lagrangian-Eulerian formulation for timedependent domains, including moving boundary layers. HMB3 uses a cell-centred finite volume approach to discrete Navier-Stokes equations on multi-block, structured grids. The 3rd order MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) approach is applied to provide high-order accuracy in space. In the present work, simulations are performed with the k-ω shear stress transport (SST) [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF]turbulence model. Regarding acoustics, the present work estimates the near-field information derived from pressure fields computed with the high-fidelity HMB3 tool. The sound pressure signal was obtained by subtracting the time-averaged pressure from the time-dependent data. All CFD grids are designed to have at least 20 cells in the near-field region to capture the target wavelength, which is calculated based on four times the Blade Passing Frequency (BPF = 400 Hz). approach has also been applied in previous work by Smith [START_REF] Smith | Acoustic analysis of counterrotating open rotors with a locked blade row[END_REF].

Summary

The numerical methodologies used in the group are summarized in Table 1. 8.

Partner

No trim is applied in all simulations. Figure 8 indicates that general tendencies captured in the measurements are reproducible using all computational tools. The differences from the test results are very much depending on the accuracy of the solvers. For example, the most of the mid-fidelity code underpredicts the torque values due to the neglection of viscosity effects. The higher order of accuracy is obtained by all CFD code except DLR TAU which underestimated the value. The choice of a vortex core model as well as a vortex core radius can have an influence on the inflow downwash velocity and therefore rotor thrust. For CFD-related methods, the proper choice of grid resolutions is key to the rotor vortex with less dissipation or dispersion errors. For this rotor, flow separations were observed in CFD results.

. Ensemble Average Velocity Fields

The simulated flow field behaviour is compared with ensemble average velocity fields obtained on 500 instantaneous velocity fields measured by the PIV technique [START_REF] Gregorio Rossignol | Multi-rotor wake interaction characterization[END_REF] as shown in Figure 10 for isolated rotor at 10120 rpm.

Figure 10 Iso-contours with streamtraces of the ensemble-average flow field colour-coded with the velocity vertical component at at Ω = 10120 RPM

The vertical velocity colour map, together with the stream traces, provides an overview of the flow field behaviour. The mean velocity shows the typical wake contraction immediately downstream of the rotor disk. To have quantitative comparison of PIV data, the comparisons are conducted on several slices at different Z positions marked as dashed line in the plot. Only vertical velocity is compared.

Figure 11 Vertical velocity trend along radial distance at Ω = 10120 RPM for different z-stations

Figure 11 shows that general forms, wake contraction and flow acceleration down to z = -1 R in the measurements are captured by most of computational tools. The comparisons are reasonably satisfactory. The Polimi DUST results exhibit possible a small dissipation model applied in the vortex development, therefore the velocity induced are very high. The results differ from other partners in large downstream position.

An interesting result is obtained in terms of the shear layer region location. Indeed, all simulations predict, more or less, the same position of the shear layer, which is different from the experiment. There is underestimation of the downward velocity in shear layer area. This means that the difference is independent of the simulation methodology used. The reason for it, is not very clear.

The difference in the area below the rotor hub occurs mainly from free wake potential code results where the hub is not modelled. Therefore, root vortex from free wake code can cause more deceleration on the flow in this area. Actually, RAMSYS and UPM results are closer to the experiment in this area. This is because the root vortex is annihilated in hover conditions. This looks a correct modelling solution. CFD results demonstrated in general similar characteristics to other results. In addition, in ONERA, UoG and ROMA3 results, as well as DLR and UoG CFD, the wind fairing Pylon is not included for the calculation, therefore large difference occurs in the root area.

Aeroacoustics

The presentation of the acoustic results will focus on data taken in the test on a polar arc represented by microphones from number 5 to 11 as defined in Figure 12. The arc is located in the Y-Z plane with an arc radius of 0.3m and polar angle difference of 15°. The microphone 9 is located in the rotational plane (X-Y), where Z=0. FW-H formulation is used by all partners.

Figure 12 Microphone positions for the isolated rotor case.

The polar directivities for the overall sound pressure level derived from microphones 5 to 11 are shown for the isolated 13x7 rotor in Figure 13. The overall sound pressure levels obtained from the test using two averaging techniques are also included. The advantage of the time-averaged noise spectrum in one rotor revolution is that the random or stochastic noise sources, such as broadband noise can be removed or reduced dramatically. For the spectrum averaging technique, the spectrum is first obtained from the time series data blocks of each rotor revolution, then averaged to obtain a mean spectrum. The comparison indicates that the acoustic simulation results are able to capture the experimental trends for all the simulation. The higher loading noise contribution for the microphone above the rotor were predicted for all partners. Very large difference occurring at microphone 5 is because the microphone five is located directly below the rotor and the rotor downwash increases the microphone self-noise, which is not simulated. In general, the contribution of the broadband noise to the OASPL is negligible in the measurements.

The noise directivity calculated by DLR, IAG and ONERA shows very good agreement. Here the comparisons involve two fidelity levels. The agreement indicates that the wake had no great influence the main noise contribution came from the loading noise, which may occur for different fidelity levels. The Polimi results closest to the experimental data, while the results ROMA3 and UoG show agreement with others. The time history data from experiment in Figure 14 low exhibits high and low frequency variations. The low frequency 2-P variations are stronger than the ones observed numerically, highlighting the discrepancies in terms of thrust already observed between the test and simulations. The computations of DLR and IAG capture well the slop when the blade approaching the microphone, but mismatch the slop when the blade pass away the microphone. The variation of the high frequency is relative weaker which could be attribute to the rotor broadband noise or some contamination due to reflections from the microphone support structure, etc.

Figure 15 compares predicted sound pressure spectra for two frequency ranges. For both the simulations and test, the rotor harmonic sound components are found to dominate at low frequency. Broadband noise becomes more important in the test for frequencies greater than 2 kHz, but the contribution to overall sound pressure level can be neglected. It should be mentioned that broadband sound and motor noise contributions were not included in numerical simulations.

Figure 15 shows all numerical simulations match very well the experimental trends, in which the highest level in the spectrum is located at the first BPF and decays almost linearly with increasing frequencies. The subharmonics among two BPFs observed in the experiment are caused by the influence of blade asymmetry and the motor noise. As all CFD simulations consider the blade asymmetry, the subharmonics are also captured.

RESULTS TWO ROTORS IN COAXIAL CONFIGURATION

For the coaxial configuration, as shown in Figure 2 and Figure 5b, the distance between two rotors is defined by ΔZ/D, where D is the diameter of the rotor in 0.33m. Both rotors are turn in opposite direction with the same rotational speed. The lower rotor defined as rotor 1, right-handed rotor (RHR) rotates counter-clockwise as for the isolated condition and the upper rotor defined as rotor 2, left-handed rotor (LHR) rotates then clockwise. The starting phase positions of the reference blade for both rotors point towards the downstream and both rotors are phase locked in the simulations. In the test however, there is no phase synchronization of the two rotors.

Ensemble Average Velocity Fields

The ensemble average velocity fields from measurement for the coaxial configuration with Δz = 0.5 R at Ω = 8025 RPM is given in Figure 17. As for the isolated case, the mean data is obtained averaging on a set of 500 instantaneous velocity map. The PIV measurement foresaw that the upper rotor was free to rotate while the lower was phase locked at rotor azimuthal position of Ψ = 90° to avoid blade laser reflection in the measurement region. The mean velocity field shows the typical wake contraction downstream the lower Rotor 1 (RHR) similar to the isolated rotor (Figure 11).

Aeroacoustics

The presentation of the acoustic results will focus on data represented by microphones on a line of array from number 0 to 5 as defined in Figure 19. The microphones are located in the X-Y plane with equal distance of 0.3m to X-Z plane. The OASPL directivity calculated by DLR, IAG, CIRA and ROMA3 shows very good agreement for all microphones although the comparisons involve two fidelity levels. Polimi results display the highest level among the simulations. All the simulations capture the trend of the directivity. Good comparisons also observed for the microphones located at positions X<=0, but large deviations from the test results for the microphones located at X>0 are observed. In order to see the contributions from individual rotors, the comparison of the individual rotor contributions from DLR simulations also included in the plot. The noise level from lower rotor 1 is in general higher than that from the upper rotor 2. The higher rotor 1 noise is caused possibly by higher non-uniform inflow introduced by the multi-inductions of the upper rotor 2 wake as well as the difference of the noise directivity. When comparing sum of two rotor noise with the individual contributions, the acoustic interferences are well visible, indicating for example clear noise cancelation for the microphones located at X>0. Due to slight asynchronization of RHR and LHR in the measurement, this phenomenon is not able to be measured.

Selected acoustic spectra obtained at two microphone positions (M3 and M4) are compared with the spectral-averaged test results, as shown in Figure 21. The spectral-averaged test data are chosen to better emphasize the relative importance of both the broadband and tonal components of the source.

Figure 21 Sound pressure level (SPL) spectra from the microphone located at 3 and 4 for 8000rpm

For the microphone located at M3, the low frequency BPF harmonics are the dominant source of noise and Rotor 1 (RHR) noise is dominant contributor to the SPL, according to DLR simulation. In the test results, there is a clear increase in tonal component levels from the 3rd BPF harmonics and sub-harmonics upward. The increase in higher harmonic tonal component levels in the test may indicate more interactions related unsteadiness especially blade leading edge interaction with as observed in [START_REF] Gregorio Rossignol | Multi-rotor wake interaction characterization[END_REF], which indicates a wider part of the blade will interact with the high-intensity part of the incoming turbulent flow caused partly by wandering of the tip vortex center. In addition, an effect of the unsynchronized speed in the experiment.

For the microphone located at M4, where large deviation from the test is observed, the test result shows BPF1 is dominant tone and large difference from the numerical simulations mainly occurs for this tone. This can be partly caused by the coherent summation used by the sum of two rotor components in the simulations. As indicated in BPF1, the sum of the contribution of the rotor 1 and 2 causes a reduction of the overall SPL due to the coherent effect. As mentioned, this phenomenon is not able to be measured due to slight asynchronization of the rotors in the measurement.

All the numerical simulations demonstrate good agreement for first 3 BPF tones. There is slightly high harmonic contribution above 2000Hz in IAG CFD results, which believe to be pure numerical artifacts.

RESULTS OF TWO ROTORS IN TANDEM CONFIGURATION

For the tandem configuration, as shown in Figure 5a, the distance between two rotors is defined by ΔZ/D and ΔX/D, where ΔZ/D=0.25 and ΔX/D=1.18 is chosen. Similar to the coaxial case, the lower rotor defined as right-handed rotor (RHR) rotates countclockwise and the upper rotor defined as left-handed rotor (LHR) rotates then clockwise.

. Ensemble Average Velocity Fields

The ensemble average velocity fields from measurement for this tandem configuration at Ω = 8025 RPM is given in Figure 23. The results were obtained on 500 instantaneous velocity fields recorded avoiding any lock between the camera acquisition frequency and the rotor speed For the tandem configuration, a better agreement in comparing with the coaxial condition between simulation and experiments is obtained. The plateau at x/R=0 is related to the modelling of the rotor hub which has been explained in previous section. The shear layer location is almost correct.

Aeroacoustics

Similar to the coaxial case, the acoustic results will be presented on a line of array from number 4 to 12 as defined in Figure 25. The microphones are located in the X-Y plane with equal distance of 0.3m to X-Z plane. The OASPL directivity pattern fit very well among all the simulations, except ROMA3 where higher noise around Rotor2 is observed. Good comparisons with the test also observed for the microphones located at positions close to rotor 1, except ROMA3 where lower noise around Rotor1 is observed. In comparison with the test, the deviation from the test results increases for the microphone towards the rotor 2. Large differences are also observed between time averaged and spectrum averaged test results.

To check the reason for the difference, selected acoustic spectrum obtained at two microphone positions (M6 and M8) are compared with the spectral-averaged test results, as shown in Figure 27.

For M6 and M8, the first 4 low frequency BPF harmonics are the dominant source of noise. For these 4 BPFs, the Rotor 1 (RHR) noise is the dominant contributor to the noise level for these selected microphones.

The large difference due to the interference effect between the simulation and test for first BPF is observed. The interference effect may not be captured as discussed before.

All the numerical simulations given in Figure 27 demonstrate good agreement for first 5 BPF tones.

Figure 1 Figure 2

 12 Figure 1 Experimental setup: single rotor (left), coaxial rotors (middle) and tandem rotors (right)

  sCMOS double-frame cameras, characterized by 5.5 Mpx and 16bit pixels dynamic range. The cameras were equipped with 50 mm Canon EF focal lenses, and mounted on motorized rotating motors, allowing for camera rotation and for Scheimplugh correction.

  3 and toward the z-axis down to z/R=-2.16 (Figure 5-b).

Figure 4 Figure 5

 45 Figure 4 PIV set-up for coaxial configuration in CIRA

Figure 7

 7 European Rotorcraft Forum, Germany, 5-7th September, 2023 This work is licensed under the Creative Commons Attribution International License (CC BY). Copyright © 2023 by author(s).

Figure 6 Figure 7

 67 Figure 6 Rotor and pylon model

Figure 8 Figure 9

 89 Figure 8 Time averaged rotor thrust and torque for 8000, 10000 and 12000 RPM

Figure 13 OASPL

 13 Figure 13 OASPL polar directivity taken from Mic 5 to Mic 11) for RPM 12000 and 10000

Figure 14 Figure 15

 1415 Figure 14 Time histories at Mic. position 9 for RPM 12000. Up: thickness noise;Low:total noise

5. 1 .

 1 Hover at RPM=8000 and Δz/D=0.255.1.1. Aerodynamic performanceA comparison of the measured and simulated time averaged rotor thrust and torque for this coaxial configuration is given in Figure16together with the value from isolated rotor. The effect of the coaxial configuration shows a decrement of the thrust. A clear effect on rotor 1, lower rotor (RHR) for both the test and all simulations (marked as red symbols) shows a larger loss of the thrust, comparing with rotor 2 LHR. This behaviour is mainly due to the effect of the upper rotor downwash which spanwise increases the induced velocity and, thus reducing the effective angles of attack.

Figure 16

 16 Figure 16 Time averaged rotor thrust and torque at coaxial configuration for 10000RPM

Figure 17

 17 Figure 17 Iso-contours with streamtraces of the ensemble-average flow field colour-coded with the velocity vertical component for the coaxial configurationwith Δz = 0.5 R at Ω = 8025 RPM

Figure 18 Figure 17 .

 1817 Figure 18 Vertical velocity trend along radial distance at Ω = 8025 RPM for different z-stations

Figure 19

 19 Figure 19 Microphone positions for the coaxial rotor case

6. 1 .Figure 22

 122 Figure 22 compares the measured and simulated time averaged rotor thrust and torque together with the value from isolated rotor. The effect of this tandem configuration shows two rotors both in numerical simulation and in the test have very close thrust and torque values. The thrust slightly lower than isolated one. All codes predict a slightly high torque value. There is no clear effect of rotor-rotor interaction observed in the averaged performance value of the coaxial configuration.

Figure 22

 22 Figure 22 Time averaged rotor thrust and torque at coaxial configuration for 10000RPM

Figure 23 Figure 24

 2324 Figure 23 Iso-contours with streamtraces of the ensemble-average flow field colour-coded with the velocity vertical component for the coaxial configurationwith Δz = 0.5 R at Ω = 8025 RPM

Figure 25

 25 Figure 25 Microphone positions for the tandeml rotor case

Figure 27

 27 Figure 27 Sound pressure level (SPL) spectra from the microphone located at 6 and 8 for 8000rpm

Table 1 Main

 1

		Code	Description
	CIRA	RAMSYS,	Unsteady, inviscid and
		ACO-suite	incompressible	free-
			wake	Boundary
			Element Method (BEM),
			Ffowcs	Williams
			/Hawkings (FW-H)
	DLR	UPM, TAU	Free	Wake	Panel
		APSIM	method, unstructured
			CFD, FW-H

Characteristics of the codes used by the partners

For a mid-fidelity aerodynamic tool based on potential formulation and free wake, the specific number of panels on the blade and pylon utilized in this paper is listed in Table

2

. The numbers are derived according to the convergence of the code results.

Table 2 Summary

 2

	Partner	Grid	Grid	Time	Num. of
		Cells	Cells	step	revs
		Rotor	Total	ISO/Multi	ISO/Multi
		(Million)	ISO/Multi
			(Million)

of the numerical resolution for mid fidelity aerodynamic tool used by the partners

The summary of CFD grids utilities is given in Table

3

.

A comparison of the measured and simulated time averaged rotor thrust and torque is given in Figure

This work is licensed under the Creative Commons Attribution International License (CC BY). Copyright © 2023 by author(s).

ACKNOWLEDGMENTS

The research leading to the presented results has been addressed within the framework of the HC/AG-26, supported by GARTEUR. POLIMI authors gratefully acknowledge the high-performance computing resources provided for this collaborative effort under the supervision of Professor Matteo Parsani by the Supercomputing Laboratory at King Abdullah University of Science & Technology (KAUST) in Thuwal, Saudi Arabia.

Copyright Statement

The authors confirm that they, and/or their company or organization, hold copyright on all of the original material included in this paper. The authors also confirm that they have obtained permission, from the copyright holder of any third party material included in this paper, to publish it as part of their paper. The authors confirm that they give permission, or have obtained permission from the copyright holder of this paper, for the publication and distribution of this paper and recorded presentations as part of the ERF proceedings or as individual offprints from the proceedings and for inclusion in a freely accessible web-based repository.