# Boundary-layer parameterization for assessing temperature and evaporation in floating photovoltaics at the utility-scale

Baptiste Amiot, Rémi Le Berre, Stéphanie Giroux-Julien, Martin Ferrand

► To cite this version:

Baptiste Amiot, Rémi Le Berre, Stéphanie Giroux-Julien, Martin Ferrand. Boundary-layer parameterization for assessing temperature and evaporation in floating photovoltaics at the utility-scale. Renewable Energy, 2026, 258, pp.124901. 10.1016/j.renene.2025.124901. hal-05405880

HAL Id: hal-05405880

https://hal.science/hal-05405880v1

Submitted on 9 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Boundary-layer parameterization for assessing temperature and evaporation in floating photovoltaics at the utility-scale

Baptiste Amiot a,*, Rémi Le Berre b, Stéphanie Giroux-Julien c,†, Martin Ferrand a,d ‡

## ARTICLE INFO

**Keywords:**

Energy meteorology

Floating photovoltaics computational fluid dynamics

Convective heat transfer coefficients

Evaporation

## ABSTRACT

A precursor model for parameterizing the effects of photovoltaic powerplants on the atmospheric boundary layer is developed using computational fluid dynamics. The method allows one to compute the surface roughness lengths, aerodynamical resistances of covered surfaces and convective heat transfer coefficients, adapted for any photovoltaic module layouts and wind directions. It has been applied for two setups: a wind tunnel system and a utility-scale floating photovoltaic installation. In these cases, the altitude-based velocity profiles was reproduced over the arrays; and we found that the turbulence generated by the photovoltaic/atmosphere interaction is greater for head- and tailwinds than sidewinds, therefore affecting the environment and the photovoltaic system. Constructing a digital twin of the floating array using large-scale meteorological fields and the parameters of the precursor model, the temperature of a monitored module was calculated and a spatial variation of 1.3 °C/km and 5.8 °C/km was estimated at the utility scale. Moreover, the waterbody evaporation was reduced by 40%-50% due to the photovoltaic panels blocking the vapour removal processes. This result decreased to 14%-20% when considering the flow spatial variations across the waterbody. Further research is necessary to adapt the parameterization to scenarios with low wind velocity.

## 1. Introduction

The energy yield from solar collectors is influenced by three factors: the surrounding meteorological conditions, the design and technology of the collectors, and the environmental interactions. The field of solar energy meteorology represents an effort to combine insights from atmospheric sciences and energy systems to improve our understanding of these elements [1], and optimize energy system outputs. Although the field relies mainly on extensive studies of radiative mechanisms, there is a growing focus on how collector alters environmental conditions, especially when considering photovoltaic (PV) systems [2]. Innovations in eco-photovoltaic systems also support this trend due to their particular interactions with the local environment and the biosphere (agrivoltaism [3], floating photovoltaic [FPV] [4]). The primary local impact of photovoltaic arrays is the shading effect, which reduces the amount of energy at the ground level, altering the temperature of the land surface [5], as well as the moisture in

the ground [6]. For FPV, shading combined with water motions and specific thermal properties reduces water temperatures [7], with additional effects on aquatic biogeochemistry [8]. A secondary impact occurs in the surroundings of the modules: the heat from the PV modules raises the air temperature [9], contributing to the heat island effect [10]. This impact of heat on the atmosphere has been found to be more related to the temperature of PV modules than to the temperature of water in the FPV context [11], highlighting the need for appropriate design of solar installations to adjust their environmental impacts [12]. In fact, the various PV layouts and the local climate modify the atmosphere/PV interaction, which ultimately adds complexity to the environmental assessment. For example in [13], a cool island was observed instead of a heating island, and it was associated with lateral airflow and local circulation patterns due to the PV array (e.g., in an environment unaltered by humans, these effects can be found in lake-desert interactions [14]). In this context, fully leveraging

<footnote>ª CEREA - Centre d'Enseignement et de Recherche en Environnement Atmospherique, Ecole nationale des ponts et chaussees, EDF RD, Marne La Vallee, 77450, France</footnote>

<footnote>ᵇ Département de Technologies et Recherches pour l'Efficacité Energetique, EDF RD, Ecuelles, 77250, France</footnote>

<footnote>ᶜ LAGEPP - Laboratoire d'Automatique, de Genie des Procedes et de Genie Pharmaceutique, Universite Claude Bernard Lyon 1, CNRS UMR5007, Villeurbanne, 69622, France</footnote>

<footnote>ᵈ Département de Mecanique des Fluides, Energie et Environnement, EDF RD, Chatou, 78401, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: baptiste.amiot@pm.me (B. Amiot).</footnote>
# Nomenclature

## Thermophysical properties

α Magnitude of the large scale variable gradients for pressure (p), temperature (θ) and scalar (y), resp. [Pa], [K], [−]

β Thermal expansion, [K⁻¹]

cₚ Specific heat capacity at constant pressure, [J K⁻¹ kg⁻¹]

Dₘ Scalar dynamic diffusion coefficient, [m² s⁻¹]

k Thermal conductivity, [W m⁻¹ K⁻¹]

p Atmospheric pressure at ground level, [Pa]

ρ Density, [kg m⁻³]

ρ₀ Density of the reference state, [kg m⁻³]

rH Relative humidity, [0–1]

T Real temperature, [°C] or [K]

θ Potential temperature, [K]

θ₀ Potential temperature of the reference state, [K]

w Specific humidity, [kg kg⁻¹]

## Kinematics & turbulence

e$^p$ Large scale variable gradients

ε Dissipation of turbulence kinetic energy, [m² s⁻³]

k Turbulence kinetic energy, [m² s⁻²]

ν Dynamic viscosity, [m² s⁻¹]

ν Turbulent viscosity, [m² s⁻¹]

τ Reynolds Stress vector, [kg m⁻¹ s⁻²]

τ₃ₖ Reynolds Stress in the k direction along a plane normal to the e₃ direction, [kg m⁻¹ s⁻²]

θₘ Wind direction, [rad] or [°]

u Velocity vector, [m s⁻¹]

u\* Friction velocity, [m s⁻¹]

uᵢ Velocity component, [m s⁻¹]

Uₘ Wind speed (norm of the velocity vector), [m s⁻¹]

z₀ Roughness length, [m⁻¹]

## Dimensionless groups & quantities

αsw Global hemispherical absorptivity (visible), [−]

εlₘ Global hemispherical emissivity (infrared), [−]

η Solar cell efficiency, [−]

FF Form factors, [−]

Prᵢ Prandtl turbulent, [−]

Π Adapted dimensionless number (ratio of the buoyancy term to the flow shear term), [−]

Y Passive scalar, [−]

## Heat fluxes & related quantities

hₐᵥ Convective heat transfer coefficient (CHTC), [W m⁻² K⁻¹]

Γ Source terms for mass m, momentum u, scalar y, energy θ), resp. [kg m⁻³ s⁻¹], [kg m⁻² s⁻²], [kg m⁻³ s⁻¹], [W m⁻³]

Lᵥ Latent heat of vaporization, [J kg⁻¹]

R Aerodynamical Resistance, [s m⁻¹]

Φ\*ₚsw Visible radiation in the PV module plane of array, [W m⁻²]

Φ\*lₘ Downward thermal radiation, [W m⁻²]

φₕ
Vertical turbulent sensible-heat flux (positive in the e₃-direction), [W m⁻²]

φₗ
Vertical turbulent flux of water vapour (positive in the e₃-direction), [kg m⁻² s⁻¹]

Uᵥ Global Transmittance, [W m⁻² K⁻¹]

## Geometry & layout

Aᵢ Surface of boundary faces i, [m²]

εₚ PV module thickness, [m]

GCR Ground coverage ratio, [m]

Hₐ
Atmosphere height, [m]

Hᵣ Floater height, [m]

Hₚ PV module height, [m]

Lᵣ Floater length, [m]

Lₚ PV module length, [m]

Sₚ PV module spacing, [m]

θₚ PV module angulation, [°]

W Numerical cell unitary width, [m]

Ωcᵢ,ᵣ Volume of numerical cell c connected to boundary face i, [m³]

## Constants

cₚ,ₐ = 1006 J K⁻¹ kg⁻¹ Specific heat capacity of dry air

cₚ,ᵥ = 1996 J K⁻¹ kg⁻¹ Specific heat capacity of dry air

g = 9.81 m s⁻² Gravitational acceleration

κ = 0.41 Karman constant

p₅ = 1 × 10⁵ Pa Reference pressure

R₉ = 287.058 J K⁻¹ kg⁻¹ Dry air gas constant

σ = 5.67 × 10⁻⁸ W m⁻² K⁻⁴ Stefan–Boltzmann constant

## Mathematical Symbols

e Unitary vector

= Ensemble average

' Fluctuation

⟨⟩ Spatial average

∂ᵢ Spatial derivative

∂ᵢ Temporal derivative

† Periodic variable tag

## Abbreviations

a Atmosphere

f Floater system

p PV module system

w Water Surface

fₚ Front surface of PV modules

rₐ Rear surface of PV modules

ABL Atmospheric Boundary Layer

CFD Computational Fluid Dynamics

CHTC Convective Heat Transfer Coefficient

PV Photovoltaic

FPV Floating Photovoltaic

C Covered Waterbody

L Land

U Uncovered Waterbody
PGD Pressure Gradient Driven

SSD Shear Stress Driven

SINF Semi-infinite Array

environmental fluid mechanics approaches are tempting, but the question of the PV system and its temperature remains.

The temperature of PV systems is usually determined from a PV thermal model and many parameterization exist in the literature [15]. They are in fact grounded to the first law of thermodynamics, assuming *steady-state* and *negligible PV module conduction*. Hence, the conceptual framework called the *big-PV module approach* reads:

$$ \alpha_{sw} \Phi_{sw}^{\text{P}} (1 - \eta(T_{sc})) = U_v (T_p - T_a), \quad (1.1) $$

where $\Phi_{sw}^{\text{P}}$ is the global irradiance in the plane of array of the PV modules, $\alpha_{sw}$ is the short wavelength absorptivity of the module, $\eta(T_{sc})$ is the solar cell efficiency related to the solar cell temperature $T_{sc}$, the so-called U-value $U_v$ is a heat loss rate that accounts for thermal radiation and convection, $T_p$ is the temperature of the module and $T_a$ the ambient air temperature. The conceptual framework is efficient in predicting module temperature for a single PV module but remains poorly adapted to PV array installations where spatial temperature gradients are observed. Such gradients are reported for various array layouts with maximum magnitudes when the array receives the most solar input (e.g., two demonstrator-sized PV racks in [16]: approximately 1 K m⁻¹ up to 10 K; single large PV racks in [17,18]: respectively 1 K m⁻¹ and 0.3 K m⁻¹, up to 10 K and 12 K). Note that spatial gradients are typically assessed over small distances and close to the corners of the PV array where the wind is highly turbulent (i.e., heterogeneous atmospheric flow region). Indeed, the wind vector contributes to the direction and magnitude of the gradient that finally influences the energy yield of the PV system [19] and the degradations in long-term performance [20]. This invites one to refine the global heat-loss rate $U_v$ using convective heat transfer coefficients (CHTC) adapted to the PV layouts [21]. However, such relations are generally not available for PV module arranged in arrays, and a supplementary effort is required to address the wind effect at scale.

The appropriate framework to simulate the PV/atmosphere interaction for utility-scale dimensions is eventually the mesoscale atmospheric model which is primarily dedicated to studying the atmospheric boundary layer (ABL). In this framework, typical physical fields of interest are the components of the wind velocity vector ($u_x, u_y, u_z$), air temperature ($T_a$) and moisture (w); their evolution and distribution are mainly governed by surface exchanges from air to ground of momentum $\tau_g$, energy $\phi_H$, and moisture $\phi_L$. Turbulence processes drive these vertical exchanges, and an appropriate scaling parameter is the friction velocity denoted $u_\star$. The proportional relation between $u_\star$ and the surface exchanges reads:

$$ \tau_g = \rho_0 u_\star^2, \quad (1.2) $$

$$ \phi_H = \rho_0 c_{p,a} u_\star \theta_\star, \quad (1.3) $$

$$ \phi_L = \rho_0 u_\star w_\star, \quad (1.4) $$

where $\rho_0$ is the air density, $c_{p,a}$ is the specific heat capacity of air, and $\theta_\star$ and $w_\star$ are scaling parameters for temperature and moisture, respectively. Fig. 1 depicts a standard ABL condition for a PV array, highlighting the strong connection between Eqs. (1.3) and (1.1). However, mesoscale frameworks are not intended to consider explicit obstacles such as buildings, trees, or PV modules; they rely on a parameterization of surface exchanges, notably for the momentum flux and the scaling parameter $u_\star$. When historically mesoscale models including PV focused on building-integrated systems for which $u_\star$ is driven by the building [22], this hypothesis does not hold when the greater element that protrudes from the surface is the PV module

itself. In order to parameterize the flow (i.e., reconstructing the vertical logarithmic evolution of wind speed above the surface), drag [23] or an equivalent roughness [24] is used and relies on the obstacle geometry. Recent studies for the configuration of PV arrays have used ad hoc roughness lengths [25], but this method requires redetermining the surface roughness from local tower measurements and is therefore bound to the layout of the PV array under investigation.

To address photovoltaic/atmosphere interactions for smaller domain sizes, computational fluid dynamics (CFD) methods are generally leveraged, as they allow fine-grained spatio-temporal resolutions. From aerodynamic studies [26], to the assessment of heat transfer for individual modules [27] or arrays [28], the size of the numerical domain remains a pain point from a computational point of view. The characteristic lengths of PV modules in front of the PV array size require a much smaller discretization length near the PVs to correctly simulate the turbulence field. A general solution is to reduce the size of the domain of interest and maintain a representative geometrical layout for the array. This is, for instance, the case in [29] where $U_v$ is computed for eighteen PV modules standing in a FPV array composed of monofloats and they determine a general relation such as $U_v = 17.7 + 5.5u$ where $u$ is the wind velocity at the inlet (in m s⁻¹). Interestingly, a temperature gradient of ≈0.4 K m⁻¹ is obtained for the first modules; however, the module temperatures stay approximately the same afterwards and the gradient exhibits a slightly negative trend (i.e., the PV module in the 18th row has a lower temperature than the one in the 6th row). Nevertheless, a new spatial discretization must be issued for any other geometrical layout, which therefore reduces the capabilities of the array slicing method. In fact, airflow patterns evolve between different layouts, directly impacting the CHTC's as found in [30] (inter-row spacings) and [31] (module heights and tilts). In addition, some characteristics of the ABL also influence the CHTC of PV modules within arrays, as shown in experimental studies under controlled conditions: the large-scale wind velocity [32] and turbulence intensity [33]. This advocates for the development of CFD methods that can handle various PV layouts without disregarding the ABL properties and their interactions with the energy yield model as described by Eq. (1.1).

This research presents a computationally tractable CFD precursor model that aims to establish a land surface parameterization of photovoltaic power plants, while providing versatility in PV design. The objective is to facilitate the assessment of operational and environmental metrics when the parameterization is integrated into a microclimate simulation (e.g., convective heat transfer coefficient, ambient temperatures, and evaporation). First, the precursor model is validated against wind-tunnel experiments. Second, a parameterization for a real industrial floating photovoltaic power plant is derived, and a microclimate model is built. Using large-scale meteorology inputs for mainly clear sky days, the temporal and spatial evolution of quantities of interest such as the temperature of PV modules and waterbody evaporation are computed, which allows us to identify the influence of the array on the environment and vice versa. Fig. 2 summarizes the workflow.

## 2. Methods

This section is segmented into three parts that introduce the proposed CFD precursor model, the microclimate CFD model supporting the PV parameterization, and the experimental survey. First, Section 2.1 is dedicated to the set of solved equations for both CFD models. Then, Section 2.2 shows the details of the precursor model and the pre-processing that allows one to determine the adapted parameters. Finally, Section 2.3 shows the CFD microclimate model using the atmospheric parameterization as well as experimental data supporting the model. In the following, Exp. 1 denotes the PV system for the geometrically resolved model mimicking the wind-tunnel experiment, whereas Exp. 2 stands for the floating PV geometry.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Illustration of the evolution of the wind-velocity profile when passing through a regularly space photovoltaic array as a result of absorption of momentum by the modules. The turbulent fluxes of energy $\phi_H$, moisture $\phi_L$ and the convective heat transfer coefficient of the PV modules, CHTC, are illustrated by thick arrows. The radiation transfers $\phi_{sw}^{net}$ and $\phi_{lw}^{net}$ are also indicated for information, they are not addressed in the work.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Illustration of the paper workflow and related sections.</figcaption></figure>

## 2.1. Generic atmospheric boundary layer model

In this article, the atmospheric model of the open source CFD solver code `saturne` [34] is used. It solves the atmospheric flow governing equations, namely the Reynolds-averaged mass, momentum, and energy balance equations, assuming incompressible flows through anelastic approximation:

$$ \partial_t (\rho_0 \bar{u}_i) = \Gamma_m, \quad (2.1) $$

$$ \partial_t (\rho_0 \bar{u}_i) + \partial_j (\rho_0 \bar{u}_j \bar{u}_j) = -\partial_i \bar{p} + \partial_j (\mu [\partial_j \bar{u}_i + \partial_i \bar{u}_j] - \rho_0 \overline{u'_i u'_j}) - \rho_0 g_i \beta (\bar{\theta} - \theta_0) + \Gamma_u, \quad (2.2) $$

$$ c_{p,d} (\partial_t (\rho_0 \bar{\theta}) + \partial_j (\rho_0 \bar{u}_j \bar{\theta})) = \partial_j (k \partial_j \bar{\theta} - c_{p,d} \rho_0 \overline{u'_j \theta'}) + \Gamma_{\theta}, \quad (2.3) $$

where $\rho_0$ denotes the mass density, $p$ is the pressure, $\underline{u}$ is the velocity vector field of components $u_x$, $u_y$ and $u_z$, and $\theta$ is the potential temperature, $g = -9.81 \epsilon_z [\text{m s}^{-2}]$ is the gravity vector, $\mu$ is the dynamic viscosity, $\beta = \frac{1}{\theta_0}$ is the thermal expansion coefficient, $k$ is the thermal conductivity of air, and $\Gamma_m$, $\Gamma_u$ and $\Gamma_\theta$ are the source terms for the mass, momentum, and energy equations, respectively. By convention, $i = 1$ (x component) refers to the East direction, $i = 2$ (y component) to the North direction, and $i = 3$ (z component) to the vertical direction.

The potential temperature is conserved over adiabatic expansion of the atmosphere with the altitude, it is defined by $\theta = T \times (\frac{\rho_s}{\rho})^{\frac{R_d}{c_{p,d}}}$, where $T$ is the temperature field, $p_s = 1 \times 10^5 \text{ Pa}$ is the reference pressure, and $R_d = 287.058 \text{ J K}^{-1} \text{ kg}^{-1}$ and $c_{p,d} = 1006 \text{ J K}^{-1} \text{ kg}^{-1}$ are respectively the gas constant and the specific heat capacity of dry air. The Reynolds-averaging operator is denoted by $\langle \cdot \rangle$ and its fluctuating part as $\langle \cdot \rangle$.

In the study, Reynolds stresses $\overline{u'_i u'_j}$ are modelled using the $k - \epsilon$ LP turbulence model [35] or the $k - \omega$ SST model [36]. A simplified gradient hypothesis is applied to model the turbulent heat flux $\overline{u'_i \theta'} = -\frac{v_i}{Pr_i} \partial_i \bar{\theta}$ in which $v_i$ is the turbulent viscosity and $Pr_i$ is the turbulent Prandtl number set to 0.7. Species mass fraction conservation equation is also solved and reads:

$$ \partial_t (\rho_0 \bar{Y}) + \partial_j (\rho_0 \bar{u}_j \bar{Y}) = \partial_j (D_m \partial_j \bar{Y} - \rho_0 (\overline{u'_j Y'}) + \Gamma_Y), \quad (2.4) $$

where $D_m$ is its molecular dynamic diffusion coefficient and $\Gamma_Y$ is a source term for Y. Note that we use the term passive scalar to describe the quantity Y in the geometrically resolved model (i.e. $Y \in [0-1]$), and the term total water content in the microclimate model. w symbol is preferred over Y symbol in this later case and $w < 1 \times 10^{-3} \text{ kg kg}^{-1}$.
<figure><img src="image_5.png" /><figcaption>Fig. 3. Sketch of the geometrically resolved CFD. Note that zone n° 4 which represents a floater in Exp. 2 vanishes in Exp. 1.</figcaption></figure>

<table><caption>**Table 1**
Dimensions of the two geometrically resolved models.</caption><thead><tr><th>Notation</th><th>Exp. 1: Validation</th><th>Exp. 2: Full scale Floating PV</th></tr></thead><tbody><tr><td>L<sub>p</sub></td><td>50.8 × 10<sup>-3</sup> m</td><td>1 m</td></tr><tr><td>B<sub>p</sub></td><td>38 × 10<sup>-3</sup> m</td><td>3 × 10<sup>-3</sup> m</td></tr><tr><td>H<sub>p</sub></td><td>63.4 × 10<sup>-3</sup> m</td><td>5.1 × 10<sup>-3</sup> m</td></tr><tr><td>ε<sub>p</sub></td><td>5.3 × 10<sup>-3</sup> m</td><td>6 × 10<sup>-3</sup> m</td></tr><tr><td>S<sub>p</sub></td><td>87.4 × 10<sup>-3</sup> m</td><td>1.44 m</td></tr><tr><td>H<sub>h</sub></td><td>8 × 10<sup>-1</sup> m</td><td>1 × 10<sup>1</sup> m</td></tr><tr><td>θ<sub>p</sub></td><td>30°</td><td>12°</td></tr><tr><td>W</td><td>5 × 10<sup>-3</sup> m</td><td>5 × 10<sup>-3</sup> m</td></tr><tr><td>L<sub>t</sub></td><td></td><td>8 × 10<sup>-1</sup> m</td></tr><tr><td>H<sub>t</sub></td><td></td><td>2 × 10<sup>-1</sup> m</td></tr></tbody></table>

Time and space discretization of the model is based on a co-located finite volume method as described in [37]. All default numerical schemes of code_saturne v8.3 are used (see [34]).

## 2.2. Geometrically resolved PV array using CFD

### 2.2.1. Computational domain and spatial discretization

The precursor model is based on a simplified computational domain that corresponds to a periodic pattern of a solar power plant, shown in Fig. 3(a). The photovoltaic module is represented by a plane of constant thickness $\varepsilon_p = 6$ mm, inclined $\theta_p$ from the horizontal and elevated $H_p$ from the ground. Other structural details are neglected due to their relatively small sizes compared to the module length. For Exp. 2, the computational domain extends over a height of 5 m and a distance $S_p = 1.44$ m in length. As shown in Fig. 3(b), portion of the fluid domain is also removed below the PV module to replicate a floating structure; the dimensions of this space are parameterized by height $H_t$ and lengths $L_t$. The spatial discretization follows a body-fitted approach with refined numerical cells close to obstacles and coarse numerical cells in the upper part of the domain. The total number of numerical cells has been designed after mesh sensitivities: 641 cells in Exp. 1 and 1763 cells in Exp. 2. Geometries and meshes are produced using the Salome platform [38] (see Fig. 4 and Table 1).

### 2.2.2. Periodic model description

The core of the methodology resides in the cyclic conditions applied to the lateral frontiers with normal directions $\varepsilon_1$ or $\varepsilon_2$, mimicking an established flow over an infinite PV array. At the top of the computational domain, the flow is driven by a shear stress vector $\tau = (\tau_{13}, \tau_{23}) = \rho_0 u_*^2 (\sin(\theta_w), \cos(\theta_w))$, where $u_*$ is the large-scale friction velocity and $\theta_w$ is the wind direction angle, which are user inputs in this section. Wall conditions are applied to the remaining obstacle surfaces,

such as the photovoltaic module, the ground, and the floater (only for Exp. 2). A wall function is prescribed for cells that are in direct contact with these surfaces. For the temperature, adiabatic conditions are prescribed at the top and bottom of the domain (including the floats in Exp. 2). A constant flux condition is set for the module surfaces $\phi_p$ to mimic the heat released from the part of the radiation that is neither converted into electricity nor remitted. Finally, Dirichlet conditions are prescribed for the passive scalar for the top and bottom of the domain and read $Y = 0$ and $Y = 1$, respectively; the module surface $Y = 0$; and for the floater surface $Y = 0$ (Exp. 2). Fig. 3(b) shows the distribution of the boundary conditions, and the applied conditions are summarized in Table 2. For the water surface, the smooth wall assumption associated with a flat surface corresponds to neglecting the effect from wave-induced motions and varying surface of exchanges (e.g., spray droplets). Hence, the methodology is limited to waterbodies with little fetch distances where waves have heights smaller than the PV system.

When the velocity field is periodic at sufficiently large downstream distance in a PV array, limiting the domain of study to a single PV arrangement requires to take into account the global pressure decrease (i.e., large-scale pressure gradient), the heat and scalar flux added to the atmosphere through the succession of PV modules (i.e., large-scale temperature and scalar gradients). To do so, we change the variables as follows:

$$p^{\#} = p - \alpha_p x_i e_i^p, \qquad (2.5)$$

$$u^{\#} = u, \qquad (2.6)$$

$$\theta^{\#} = \theta - \alpha_0 x_i e_i^p, \qquad (2.7)$$

$$Y^{\#} = Y - \alpha_Y x_i e_i^p, \qquad (2.8)$$

where $\alpha_0^{\#}$ denotes the periodic variables, $\alpha_0$, $\alpha_Y$ and $\alpha_p$ are the magnitude of the large-scale gradients in the direction $e^p$. Injecting Eqs. (2.5)–(2.7) in Eqs. (2.1)–(2.3) gives the following system of equations that verifies periodic conditions:

$$\partial_t (\rho_0 \bar{u}_i^\#) = 0, \qquad (2.9)$$

$$\begin{aligned} \partial_t (\rho_0 \bar{u}_i^\#) + \partial_j (\rho_0 \bar{u}_j^\# \bar{u}_j^\#) = & - \partial_l \bar{p}^\# \\ & + \partial_j (\mu [\partial_j \bar{u}_i^\# + \partial_l \bar{u}_j^\#] - \rho_0 \overline{u_j^\# u_l^\#}) \\ & - \rho_0 g_s \beta (\bar{\theta}^\# - \theta_0) \underbrace{- \alpha_p e_i^p}_{\Gamma_{u^\# , i}} \end{aligned} \qquad (2.10)$$

$$c_{p,d} \left[ \partial_t (\rho_0 \bar{\theta}^\#) + \partial_j (\rho_0 \bar{u}_j^\# \bar{\theta}^\#) \right] = \partial_j \left( k \partial_j \bar{\theta}^\# - c_{p,d} \rho_0 \overline{u_j^\# \theta^\#} \right)$$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Spatial discretization of the two geometrically resolved models.</figcaption></figure>

<caption>**Table 2**</caption>

Summary of boundary conditions applied to the total volume of control. The number of the boundaries are indicated in Fig. 3(b). Note that the conditions exerted by smooth walls are implemented with wall functions (two scales smooth-rough wall functions).

<table><thead><tr><th>Surface name</th><th>n°</th><th>Condition</th><th>Velocity</th><th>Temperature</th><th>Passive scalar</th></tr></thead><tbody><tr><td>Fictitious surface</td><td>1<br>1'</td><td>Periodic</td><td>u<sub>1</sub>|<sub>x=0</sub> = u<sub>1</sub>|<sub>x=s<sub>p</sub></sub><br>u<sub>1</sub>|<sub>x=s<sub>p</sub></sub> = u<sub>1</sub>|<sub>x=w</sub></td><td>θ<sub>1</sub>|<sub>x=0</sub> = θ<sub>1</sub>|<sub>x=s<sub>p</sub></sub><br>θ<sub>1</sub>|<sub>x=s<sub>p</sub></sub> = θ<sub>1</sub>|<sub>x=w</sub></td><td>Y<sub>1</sub>|<sub>x=0</sub> = Y<sub>1</sub>|<sub>x=s<sub>p</sub></sub><br>Y<sub>1</sub>|<sub>x=s<sub>p</sub></sub> = Y<sub>1</sub>|<sub>x=w</sub></td></tr><tr><td>Water</td><td>2</td><td>Smooth wall</td><td>u = 0 m s<sup>-1</sup></td><td>φ = 0 W m<sup>-2</sup></td><td>Y = 1</td></tr><tr><td>PV Module</td><td>3</td><td>Smooth wall</td><td>u = 0 m s<sup>-1</sup></td><td>φ<sub>F</sub> = 200 W m<sup>-2</sup></td><td>φ<sub>Y</sub> = 0</td></tr><tr><td>Float</td><td>4</td><td>Smooth wall</td><td>u = 0 m s<sup>-1</sup></td><td>φ = 0 W m<sup>-2</sup></td><td>φ<sub>Y</sub> = 0</td></tr><tr><td>Atmosphere</td><td>5</td><td>Inlet</td><td>τ<sub>13</sub> = sin (θ<sub>w</sub>) ρ<sub>0</sub>u<sub>1</sub><sup>*</sup><br>τ<sub>23</sub> = cos (θ<sub>w</sub>) ρ<sub>0</sub>u<sub>1</sub><sup>*</sup></td><td>φ = 0 W m<sup>-2</sup></td><td>Y = 0</td></tr></tbody></table>

$$ \underbrace{-c_{p,d} \rho_0 \alpha_{\theta} \bar{u}_j^\sharp e_j^p}_{\Gamma_{\theta^\sharp}} \quad (2.11) $$

where the source terms $\Gamma_{u^\sharp,j}$ and $\Gamma_{\theta^\sharp}$ depend on the type of boundary condition prescribed. For shear-driven flows, the contribution from $\alpha_p$ in momentum equation cancels out ($\Gamma_{u^\sharp,j} = 0$). Note that the periodic contribution of the buoyant term in the momentum equation is balanced by the hydrostatic part of the pressure field. The large-scale variable gradient for the potential temperature $\alpha_\theta$ is computed from the heat balance for a single pitch, as heat flux at the top of the domain is set to zero (see Table 2):

$$ \alpha_{\theta} = \frac{2 \times \phi_{p} (\mathcal{L}_{p} + \varepsilon_{p})}{c_{p,d} \rho_{0} \int_{\Omega_{tot}} \bar{u}_{i}^{\sharp} e_{i}^{p} d\Omega} \quad (2.12) $$

where the volume of the computational is $\Omega_{tot}$. Applying mass balance on $Y^\sharp$ with boundary conditions considered in Table 2 gives $\alpha_Y = 0$ and $\Gamma_Y = 0$ in Eq. (2.4).

Exp. 1 requires approximately 160 s of physical time to reach the steady state when initialized with a log-law velocity profile, a uniform temperature of 15 °C, and a uniform scalar value of 0. In Exp. 2, physical time to reach the steady state varies from 600 s to 1500 s, depending on velocity frictions, wind directions, and the height of the computational domain. Small wind frictions, transversal winds, or larger domain

heights increase this time. In this configuration, a domain size of 5 m is used because it allows us to replicate the flow structures obtained for larger domain sizes and reduces the time to solution.

### 2.2.3. Momentum, moisture and energy parameterizations

Periodic airflow quantities exhibit a logarithmic profile in the upper part region of the domain (z in $[\gamma, H_a]$) when steady state is reached, and we use this domain to determine a parameterization of the vertical turbulent fluxes Eqs. (1.2), (1.3) and (1.4). For clarity, we simplify the notation by omitting the ensemble average and periodic symbols in this subsection.

First, parameterizing the momentum absorption by the PV array requires defining the average profile of large scale velocity as:

$$ \langle u_{\perp} \rangle (z) \approx \frac{u_{\star}}{\kappa} \ln \left( \frac{z}{z_0} \right), \quad (2.13) $$

where $u_\perp = u_x \cos(\theta_w) + u_y \sin(\theta_w)$ is the velocity in the large scale wind direction $\theta_w$ and $z_0$ is the roughness length to be determined. The symbol $\langle \cdot \rangle$ denotes a spatial average for a given slice of atmosphere at height z which reads:

$$ \langle u_{\perp} \rangle (z) = \frac{1}{S_p \times W} \int_0^{S_p} \int_0^W u_{\perp} dy dx, \quad (2.14) $$

where W is the cell width (set to 1 m). As the friction velocity is prescribed by the top boundary condition, the roughness length $z_0$ of
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Sketch of parameterization for intensive scalars where Y is either the potential temperature field θ or the humidity scalar Y. We used γ = 2H_p for θ and γ = 0.5H_p for Y.</figcaption></figure>

the average velocity profile can be obtained by minimizing the mean square error in the upper domain region (γ = 2H_p), given by:

$$\epsilon(z_0) = \int_{\gamma}^{H_a} \left( \langle u_{\perp} \rangle(z) - \frac{u_{\star}}{\kappa} \ln \left( \frac{z}{z_0} \right) \right)^2 dz, \quad (2.15)$$

which finally leads to:

$$\ln(z_0) = \frac{\sum_{c/z_c > \gamma} \left( \ln(z_c) - \kappa \frac{u_{\perp c}}{u_{\star}} \right) \Omega_c}{\sum_{c/z_c > \gamma} \Omega_c}, \quad (2.16)$$

where Ω_c is the volume of the numerical cell c. Then for the moisture parameterization, Eq. (1.4) is rewritten using the aerodynamic resistance R_a between the bottom of the domain and the spatial average value 〈Y〉_up, as illustrated in Fig. 5(a), such that:

$$\phi_L = - \frac{\rho_0}{R_w} (Y_w - \langle Y \rangle_{up}). \quad (2.17)$$

〈Y〉_up is defined by:

$$\langle Y \rangle_{up} = \frac{1}{S_p \times W \times (H_a - \gamma)} \int_{\gamma}^{H_a} \int_{0}^{S_p} \int_{0}^{W} Y dxdydz. \quad (2.18)$$

As Y_w is prescribed by the boundary conditions, φ_L and 〈Y〉_up are deduced from the computation which gives R_a.

Finally, the CHTC for the front surface of the PV module is deduced from the computation of T_p by the fluid solver, see Fig. 5(b). It reads:

$$h_{cv,fr} = \frac{\phi_p}{(T_p - \langle T \rangle_{up})}. \quad (2.19)$$

## 2.3. Microclimate model and observational survey

### 2.3.1. Digitalization of the environment

The atmospheric parameterization is integrated within a microclimate model designed to simulate the environment of an industrial FPV powerplant located in the Hautes-Alpes region of southern France (44.34N, 5.86E).

Fig. 6 shows the power plant occupying two-third of a waterbody, as well as the nearby environment dominated by crop and forest lands with distant high-mountain terrain. In this regard, the simulation domain has been constructed using a cylinder of 2.35 km radius whose bottom surface matches the topography, see Fig. 7(a). Buildings, trees and crops are incorporated as surface roughness elements based on the land-use calculated from the satellite imagery (Fig. 6(b)), and the most likely values for the roughness lengths are applied, see [39]. The

space occupied by the photovoltaic array is reproduced from industrial sketches rather than satellite imagery to maximize geometric precision, and a dedicated roughness length is assigned (see Section 2.3.3). Figs. 7(b) and 7(c) show the result for the western island of the powerplant. Then, the height of the domain is set to 1200 m and the numerical discretization is refined in the lower region of the domain (i.e., approximately 4 m² per 5 m height around the PV islands). Considering larger numerical cells in the upper part of the domain and around the lateral sides (i.e., approximately 80 m² per 50 m height), the mesh totals 1.1 × 10⁶ cells approximately.

### 2.3.2. Characteristics of the computation and environmental dataset

Eqs. (2.1)–(2.4) hold in the numerical domain assuming that the thermodynamic variables are the total water content Y → w and the liquid potential temperature θ → θ_l. The total water content w represents the sum of liquid water w_l and vapour w_v, with phase transitions modelled through subscale models. For the liquid potential temperature θ_l, the approximation proposed in [40] is applied: θ_l = θ - L_v/θ_c w_l, since the simulated meteorological conditions ensure that L_v/w_l/c_p dT ≪ 1 [41]. This allows to remain conservative with respect to the thermodynamic processes, as well as simplifying the numerical resolution. Note that the source terms in all the equations are set to zero for all cells in the domain, with the exception of the cells in direct contact with the solar array surfaces; this is further developed in Section 2.3.3.

The microclimate model is then run for various large-scale meteorological conditions based on hourly meteorological records provided by the Météo-France agency at the B1 site (wind velocity and direction at 10 m height, air temperature, and relative humidity). The environmental quantities and PV module temperature recorded at A1 and M1 are used for validation purposes (see Table 3). At the A1 location, a floating raft embedding meteorological probes was specifically designed and mechanically connected to the floating powerplant, with supplementary mooring lines and gravity anchors to minimize motion. The probe set is composed of an ultrasonic anemometer (10 Hz), a weather station (air temperature T_a, relative humidity rH_a, pressure p_a), and a long-wave radiative sensor for downward flux φ_{lw}⊥. The anemometer is positioned to reduce the sheltering effects of the other materials and probes (Fig. 8(a)); however, it is assumed that only the east winds up to the south winds are representative of the airflow in the array. Winds recorded from the west and north directions are in fact characteristic of the nearby crop lands and the mineral embankment (2.5 m height approximately). Details on the anemometer dataset and treatments are discussed in Appendix A.2. Hence for the sake of representativeness, only easterly winds are considered in this study. At the M1 location,
<figure><img src="image_8.png" /></figure>

(a) Sensor set locations including Meteo-France weather station (B1). Satellite pictures from Copernicus and USGS.

<figure><img src="image_9.png" /></figure>

(b) Land-use distribution retrieved from THEIA database (10 m² resolution) and based on classes defined in [39].

<caption>Fig. 6. Satellite views of the microclimate study case.</caption>

<caption>**Table 3**</caption>

Details of installed sensors on the FPV powerplant as well as the data collected by the Meteo-France Weather station.

<table><thead><tr><th>Location</th><th>Instrument</th><th>Parameter</th><th>Sampling timestep (s)</th></tr></thead><tbody><tr><td rowspan="4">A1</td><td>SGR4</td><td>Φ<sub>lw</sub><sup>†</sup></td><td>60</td></tr><tr><td>GMX600</td><td>T<sub>a</sub>, rH<sub>a</sub>, p<sub>a</sub></td><td>60</td></tr><tr><td>PT-100</td><td>T<sub>v</sub></td><td>60</td></tr><tr><td>3-D Anemometer</td><td>u<sub>x</sub>, u<sub>y</sub>, u<sub>z</sub></td><td>0.1</td></tr><tr><td rowspan="2">M1</td><td>MS80S</td><td>Φ<sub>lw</sub><sup>‡p</sup></td><td>60</td></tr><tr><td>PT-100</td><td>T<sub>p</sub></td><td>60</td></tr><tr><td>B1</td><td>Meteo-Station</td><td>U<sub>w</sub>, rH, T, θ<sub>w</sub></td><td>3600</td></tr></tbody></table>

a short-wavelength radiation sensor in the plane-of-array and five PT-100 sensors affixed to the rear side of two PV module are also installed (Figs. 8(b) and 8(c)). Combined with the downward thermal radiation, the recorded plane-of-array irradiation directly feeds the PV model into the microclimate model, so that the uncertainty in the computation of the PV module temperature is principally due to the wind effect and the CHTC. The measurements from the five temperature probes are averaged out and therefore used to assess model accuracy.

This leads to considering 72 hourly situations that correspond to the number of easterly wind occurrences, observed at A1, during a two-month period in January and February 2025, and for which the plane-of-array irradiance, measured at M1, is greater than 250 W m⁻². A specific day of interest in August 2023 is also used because it features a large number of successive easterly wind occurrences with clear-sky conditions typical of a summer period, see Fig. 9. Therefore, the design of experiment consists of 72 independent steady-state simulations using a 1 s time step during 6000 s per simulation. Note that for each simulation, all computed fields (including module temperature) are temporally averaged between iterations 4000 s and 6000 s.

### 2.3.3. Boundary conditions and treatment of waterbody surfaces

Each hourly simulation uses the records at the B1 location to determine large-scale meteorological fields (i.e., velocity, temperature, and humidity), which are then applied as boundary conditions using neutral ABL assumptions. Importantly, the inlet or outlet nature of the boundary conditions prescribed to the vertical sides of the domain is based on the flow direction. For inlets, altitude-based Dirichlet conditions are set:

$$u(z) = \frac{u_{\star}}{\kappa} \ln \left( \frac{z}{z_0} \right) [-\sin(\theta_w)\epsilon_1 - \cos(\theta_w)\epsilon_2], \quad (2.20)$$

with $z_0 = 0.01$ m and $u_{\star}$ recomputed from the wind velocity records. Turbulence profiles are set using $k = \frac{u_{\star}^2}{\sqrt{C_{\mu}}} (C_{\mu} = 0.09)$ and $\epsilon = \frac{u_{\star}^3}{\kappa(z+z_0)}$.

For outlets, homogeneous Neumann conditions are imposed for all variables, except for the pressure which is imposed. When the flow becomes parallel to an outlet surface, homogeneous Neumann conditions shift to Dirichlet-type boundaries to maintain numerical stability. Note that a symmetry condition is prescribed at the top of the domain.

Then, the set of boundary conditions applied to the bottom of the domain is differentiated between land (L), free water surface (FS) and surface covered by PV (PV). Considering the parameterization on momentum, the three surfaces prescribe a two-scale wall function based on roughness lengths, respectively determined from: land cover categorization (L), the Charnock law from [42] (FS):

$$z_0 = 1.2 \times 10^{-2} \times \frac{u_{\star}^2}{g}, \quad (2.21)$$

and the wind-dependent relation elaborated in Eq. (2.16) (PV).

For temperature boundary conditions, the land (L) and the covered water (PV) are associated to no-flux conditions whereas the free water (FS) is set to a Dirichlet condition using the measured water temperature $T_w$. Finally, for total water contents, free water (FS) and covered water (PV) are associated to Dirichlet-type and Robin-type boundary conditions, respectively. In both cases, the Dirichlet contribution is applied using the empirical correlation from [43]:

$$w(T_w) = 1.35 \times 10^5 \times \exp \left( \frac{-4157}{T_w - 33.91} \right), \quad (2.22)$$

where $T_w$ is expressed in K. The Robin-type boundary condition is prescribed by forcing the resistance $R_w$ computed in Eq. (2.17). A no-flux condition is also prescribed for land surfaces (L).

The numerical cells in direct contact with the covered and free waterbody must take into account the effect of evaporative transfers because the water temperature is set (FS) or assumed (PV) constant. This contribution is implemented using a source term in energy conservation Eq. (2.3). Formally, let $A_i$ be the surface of a face i and $\Omega_{c,i}$ be the volume of the cell c adjacent to the surface i, the source term $\Gamma_{\theta_i}$ to implement is composed of three contributions:

$$\Gamma_{\theta_i} = \underbrace{(\Gamma^{\text{sens}} + \Gamma^{\text{lat}})}_{\text{From evaporation transfer (FS+PV)}} + \underbrace{\Gamma^{\text{P}}}_{\text{From module release (PV)}}, \quad (2.23)$$

where $\Gamma^{\text{sens}}$ is the sensible contribution from evaporated water, $\Gamma^{\text{lat}}$ is the latent heat contribution, and $\Gamma^{\text{P}}$ is the sensible heat released by the module for only the covered faces (PV). First, the sensible heat from the evaporation process is determined considering that the evaporated mass $w_i - w_{c,i}$ gives an energy to the atmosphere equivalent to the difference in enthalpy between $T_i$ and $T_{c,i}$. Assuming that the
<figure><img src="image_10.png" /><figcaption>(a) Spatial discretization of the computational domain, cell sizes are finer at the domain centre where the PV array is located, and coarser towards the domain edges.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Satellite view of the floating photovoltaic (c) Horizontal cross-section of the mesh at water powerplant (west side).</figcaption></figure>

<caption>Fig. 7. Spatial discretization retained for the microclimate case study and illustration of the digitization for one floating photovoltaic island of the array as delineated face areas in the spatial discretization.</caption>

specific heat capacity of the vapour $c_{p,v}$ is weakly dependent on the temperature, the given energy becomes $c_{p,v} \times (T_i - T_{i,c})$. Therefore, the first contribution reads:

$$ \Gamma^{\text{sens}} + \Gamma^{\text{lat}} = \left(\frac{\rho_0}{R_u}\right) (w_i - w_{c,i}) [c_{p,v}(T_i - T_{c,i}) + L_v] \frac{A_i}{\Omega_{c,i}}, \quad (2.24) $$

where $R_u$ is the aerodynamical resistance: of the covered water described in Eq. (2.17) for the covered faces (PV), of the open-water using the Charnock roughness length described in Eq. (2.21) for the free surface (FS), and $L_v$ is the latent heat of vaporization.

Finally, the term $\Gamma^P$ in Eq. (2.23) is only applied for covered faces (PV) and it requires to describe the evolution of the PV module temperature $T_p$. We do not use the conceptual framework described in Eq. (1.1) but an alternative version that uses the computed CHTC (Eq. (2.19)). A description of the model is given in Appendix A.1. The heat

transferred through convective mode then reads:

$$ \Gamma^P = \text{GCR}(h_{cv,fr} + h_{cv,\text{re}}) (T_p - T_{c,i}) \frac{A_i}{\Omega_{c,i}}, \quad (2.25) $$

where $\text{GCR} = \frac{L_p}{S_p}$ is the ground coverage ratio, and the CHTC for the rear face of the PV module is set to $h_{cv,\text{re}} = 1 \text{ W m}^{-2} \text{ K}^{-1}$.

## 3. Results

### 3.1. Flow dynamics of the geometrically resolved model

First, the precursor model is applied to Exp. 1 design in order to validate the proposed numerical approach. Two types of flow forcing for the precursor model are examined: shear-stress driven (SSD, $u_\star$)
<figure><img src="image_12.png" /><figcaption>Fig. 8. Pictures of the main instruments for the microclimate field campaign. General location is indicated in Fig. 6(a).</figcaption></figure>

applied to the upper boundary condition $\alpha_p = 0$) and pressure gradient driven (PGD, $\alpha_p \neq 0$ and the top boundary becomes a smooth-rough wall condition). Furthermore, a semi-infinite array set-up is constructed by replicating Exp. 1 a hundred times in the streamwise direction, it is labelled SINF. Fig. 10 shows that the SSD, PGD, and SINF models successfully represent the velocity distributions of the wind tunnel experiment. Among these, the SINF model exhibits the minimal error margin and captures the velocity inflection $z/H_p \approx 0.95$, suggesting that the RANS approach and the $k-\omega$ SST turbulence model are adapted to simulate photovoltaic arrays interacting with the atmosphere. However analysing SINF, one can see that the velocity profile has not reached an invariant state by pitch translation at the 9th row especially for the subpanel flow: this state is reached around the 33th rows matching the PGD velocity profile. Larger recirculation areas and lower magnitude in velocities around the rear side of the obstacle are observed when permanent conditions are reached, translating into a lower altitude for the inflection $z/H_p \approx 0.6$. From these results, one can notice that the velocity profile under the panel depends on the way the flow is driven (pressure-driven or shear driven) and the match with wind tunnel experiment is not good in this area, which is probably due to the fact that the flow is not periodic at this range.

Although the method to drive the atmosphere flow influences the wind velocity below the PV module, it could be observed that the velocity inflection heights (altitude where $u_y = 0$) are identical between SSD and PGD approaches and the wind profiles above the PV panels are similar. That suggests that momentum exchanges between the air above the panel and below are well reproduced, and so should be the CHTC and the aerodynamical resistance for vapour. As SSD is coherent with the atmospheric surface layer assumption (constant fluxes), this numerical approach is preferred for the geometrically resolved model.

## 3.2. Parameterization of floating photovoltaics

To assess the influence of the FPV configuration on atmospheric quantities (Exp. 2), we calculate the dynamical roughness length, the aerodynamic resistance of the covered water and the CHTC of the

PV panel. Seven simulations are performed with an evolving wind direction and fixed shear stress set to $u_* = 0.36 \text{ m s}^{-1}$. Note that for the fully transverse wind ($\theta_w = 90^\circ$) a 3-D precursor simulation was used. Fig. 11(a) highlights that tailwinds and headwinds yield large roughness lengths, while sidewinds minimize them; these indicate that the geometry of photovoltaic arrays plays a significant role in influencing momentum absorption. Substantial evidence supporting the model's accuracy is demonstrated by the headwind results, where numerical and empirical results are aligned. The evolution in the magnitude of the roughness is also notably consistent: computational estimates range from $1 \times 10^{-1} \text{ m}$ and $2 \times 10^{-3} \text{ m}$, and experimental data exhibit a similar with a greater roughness length for headwinds, whereas side winds are associated with a lower roughness length. An explanation for the larger magnitudes of the roughness lengths for side winds in the experiment is the fluctuating wind direction: each roughness length measurement considers a range of wind directions with distinct PV-atmosphere interactions. In the precursor model, the wind direction fluctuation does not exist and, therefore, the wind-surface interaction is minimal. To ease future utilization of the numerical results, the roughness length is fitted using the following form:

$$ \left\{ \begin{array}{ll} z_0 = a_1 \theta_w^2 + b_1 \theta_w + c_1, & \text{for } \theta_w \in [0; \pi[, \\ z_0 = a_2 \theta_w^2 + b_2 \theta_w + c_2, & \text{for } \theta_w \in [\pi; 2\pi[, \\ \text{but } z_0 = z_0(\theta_w = \pi/2), & \text{if } z_0 \le z_0(\theta_w = \pi/2), \end{array} \right. \quad (3.1) $$

where the set of coefficients $a_1, \dots, c_2$ reads:

$$ \begin{bmatrix} a_1 & b_1 & c_1 & a_2 & b_2 & c_2 \\ 0.03457 & -0.11675 & 0.09268 & 0.03457 & -0.31478 & 0.70536 \end{bmatrix} = \begin{bmatrix} 0.03457 & -0.11675 & 0.09268 & 0.03457 & -0.31478 & 0.70536 \end{bmatrix} \quad (3.2) $$

and $z_0(\theta_w = \pi/2) = 0.002 \text{ m}$.

A second element that can affect the roughness length is the buoyancy effect due to PV heating. We evaluated the ratio of turbulence due to buoyancy relative to that due to shear using a Richardson-like number defined by:

$$ \Pi = \frac{\beta g \phi_p L_p}{\rho_0 c_{p,d} u_*^2}. \quad (3.3) $$
<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Environmental conditions for the selected observation day. The dataset is averaged on hourly time steps and data are associated to the beginning of each time interval. It includes measurements from observation points above the lake (A1, A3, M1) and at the Meteo-France weather station located at 5 km from the lake, denoted B1. Note that the A3 point of measure is located on the far east of the lake and is deemed to be not perturbed by the PV array.</figcaption></figure>

It is observed that the flow becomes laminar when $Π > 1$ and remains fully turbulent otherwise; therefore, we tested this hypothesis by simulating new situations with height prescribed heat fluxes and keeping shear stress constant ($u_* = 0.36 \text{ m s}^{-1}$). The results are presented in Fig. 11(b) and yield an exponential relation $z_0|_{\theta_w=0^\circ} = z_0(0) \times \exp(0.7\Pi)$ with $z_0(0) = 1.1 \times 10^{-1} \text{ m}$. Under normal conditions ($\Phi_{sw}^{\lambda p} = 1000 \text{ W m}^{-2}$, $T_a = 20^\circ\text{C}$, $U_{\parallel}|_{z=3 \text{ m}} \in [2; 5] \text{ m s}^{-1}$), the Richardson-like number lies between $1 \times 10^{-1}$ and $1$ and the dynamical roughness length is increased between 10% and 100%. In consequence, the roughness length obtained by the CFD precursor model may underestimate the real roughness length, which is in line with the observed trend especially for side winds. Therefore, the exponential relation limits the precursor model when the shear is small ($Π \gg 1$) and the numerical relation for the roughness length (Eq. (3.1)) may have discrepancies.

Finally, the aerodynamic resistance of the covered water and the CHTC PV panel is presented in Fig. 12. The direction of the wind significantly modifies the transfer rates. For evaporation, the mass transfer rate is greater for tailwind conditions, which is explained by

the channelling effect below the photovoltaic module, which allows the atmosphere to stir vapour efficiently compared to the headwind case where the photovoltaic panel blocks the transfer and acts as a mechanical buffer. For CHTC, there are no such differences with respect to the wind direction, and both tail and head winds have a comparable influence on the heat transfer rates; only the side winds are prone to lower heat exchanges, potentially because of the little mechanical interaction in these circumstances. With the exception of aerodynamic resistance under tailwind flows, the transfer rates show a similar evolution to the roughness length. Therefore, the highest evaporative rate would be fed by a more turbulent flow, increasing the global evaporative transfer. The stability of the atmosphere could also influence the transfer rates as shown in Fig. 11 for momentum. As warm air is lifted upward due to natural convection, the aerodynamic resistances (resp. CHTC) would decrease (resp. increase), with varying magnitudes based on the rate of buoyancy over mechanical shear. This could be taken into account using the precursor model by adding one degree of freedom in the design of simulations (e.g., wind directions
<figure><img src="image_14.png" /><figcaption>Fig. 10. Velocity profiles for the Exp. 1 case, compared to the experimental results assuming a reference velocity of $U_0 = 3.5 \text{ m s}^{-1}$.</figcaption></figure>

<figure><img src="image_15.png" /></figure>

(a) Directional roughness lengths from the precursor model and the experiment.

(b) Sensitivity of the roughness length to the buoyancy ($\theta_w = 0^{\circ}$, N).

<caption>Fig. 11. Evaluation of the roughness length for the FPV design (Exp. 2). Dashed lines represents the best fit for the numerical results. Boxplots give an indication on the dynamical roughness experiences at the industrial FPV site.</caption>

and stability functions, as the proposed Eq. (3.3)), thus increasing computational effort. Opting for purely mechanically-driven turbulence is a strategic decision to integrate simple transfer rate functions into the microclimate model. The ones adapted for the FPV system reads:

$$ \begin{cases} y = a_1 \theta_w^4 + b_1 \theta_w^3 + c_1 \theta_w^2 + d_1 \theta_w + e_1, & \text{for } \theta_w \in [0; \pi[, \\ y = a_2 \theta_w^4 + b_2 \theta_w^3 + c_2 \theta_w^2 + d_2 \theta_w + e_2, & \text{for } \theta_w \in [\pi; 2\pi[, \end{cases} \quad (3.4) $$

$$ A_{h_{cv,fr}} = \begin{bmatrix} a_1 & a_2 \\ b_1 & b_2 \\ c_1 & c_2 \\ d_1 & d_2 \\ e_1 & e_2 \end{bmatrix} = \begin{bmatrix} 0.3 & 0.3 \\ -2.0 & -5.5 \\ 3.5 & 36.0 \\ -0.4 & -101.0 \\ 7.3 & 108.0 \end{bmatrix} \quad (3.5) $$

and:

$$ A_{R_y} = \begin{bmatrix} a_1 & a_2 \\ b_1 & b_2 \\ c_1 & c_2 \\ d_1 & d_2 \\ e_1 & e_2 \end{bmatrix} = \begin{bmatrix} -9.1 & -9.1 \\ 65.5 & 162.0 \\ -146.2 & -1055.5 \\ 108.1 & 2951.8 \\ 48.4 & -2902.3 \end{bmatrix} \quad (3.6) $$

Treating water content as a transport variable in the precursor model could also extend its applicability, particularly in edge cases where surface condensation and evaporation take place on the PV module (e.g., early in the morning). Under such conditions, both humidity
<figure><img src="image_16.png" /><figcaption>Fig. 12. Evaluation of the aerodynamic resistances of the covered water and the convective heat transfer coefficient for the front surface of the PV module, as a function of the wind direction (velocity friction set constant: $u_{\infty} = 0.36 \text{ m s}^{-1}$).</figcaption></figure>

<table><caption>**Table 4**
Configuration of simulation study cases and numerical performances for the 72 h of eastern winds.</caption><thead><tr><th rowspan="2">Label</th><th rowspan="2">Thermal model (PV)</th><th rowspan="2">Ambient Temp.</th><th rowspan="2">Wind vector</th><th rowspan="2">CHTC</th><th colspan="2">T<sub>p</sub> (M1)</th><th colspan="2">&part;<sub>i</sub>T<sub>p</sub></th></tr><tr><th>MBE</th><th>RMSE</th><th>MEAN</th><th>MAX</th></tr></thead><tbody><tr><td>Baseline 1</td><td>Standalone</td><td>Measured</td><td>Measured</td><td>literature</td><td>-6.1 &deg;C</td><td>6.8 &deg;C</td><td>-</td><td>-</td></tr><tr><td>Baseline 2</td><td>Standalone</td><td>Measured</td><td>Measured</td><td>Optimized</td><td>0 &deg;C</td><td>2.4 &deg;C</td><td>-</td><td>-</td></tr><tr><td>Case 3</td><td>Standalone</td><td>Measured</td><td>Measured</td><td>Computed</td><td>0.2 &deg;C</td><td>2.6 &deg;C</td><td>-</td><td>-</td></tr><tr><td>Case 4</td><td>In CFD</td><td>Computed</td><td>Computed</td><td>Computed</td><td>0.5 &deg;C</td><td>3.0 &deg;C</td><td>1.3 &deg;C</td><td>5.8 &deg;C</td></tr></tbody></table>

and temperature fields influence energy conservation, requiring a supplementary degree of freedom in the design of simulations. Although buoyancy-driven changes would also modify transfer rates, this impact is likely minor compared to the temperature gradient between the PV system and the surrounding air. Consequently, for classical conditions in FPV operation, omitting this effect remains justified.

### 3.3. Temperature of floating photovoltaics at the array scale using numerically obtained CHTC relation

The boundary-layer parameterization is then implemented into the microclimate and PV models to simulate the temperature of the PV module at M1, and the spatial temperature gradient in the array. Four modelling configurations with varying data inputs and parameter choices are established. The first three scenarios only use the PV model described in Appendix A.1, with measured ambient temperature and wind speeds as inputs, but considering distinct CHTC relations: characteristic of well-cooled PV system $h_{cv,fr} = 8.56 + 2.55U_w$ as observed in [44], denoted Baseline 1; optimized from the experimental dataset, $h_{cv,fr} = 1.5 + 1.5U_w$ and denoted Baseline 2; from the ABL parameterization, Eqs. (3.4) and (3.5) and labelled Case 3. The last scenario denoted Case 4 is the microclimate model where the atmosphere and the PV models are solved using the full ABL parameterization.

Table 4 summarizes the results for the four configurations. First comparing Baseline 1 with Baseline 2 demonstrates that the CHTC is essential for accurate FPV temperature modelling and that well-cooled CHTC relationship is not applicable to mono-float FPV layouts. These results extend the observation made in [45] that wind exposure is the critical parameter in the cooling effect of FPV systems, for utility-sized FPV installation. As no effect related to the waterbody environment is modelled (e.g., evaporative cooling or radiative cooling from the cooler water surface), it is hypothesized that these effects would remain marginal in the cooling of mono-float FPV system; the geometry

of the PV system is likely to be the main criteria to play with in order to achieve notable cooling effect for FPVs. As shown in Case 3 and the small bias of 0.2 °C, the biperiodic model appears to be an interesting solution to determine this relation between the CHTC and the PV module/array layouts (e.g., module height, tilt and pitch). Because geometry can be easily altered, this approach can be used to evaluate different PV geometries and their heat transfer performance, especially during the energy resource assessment phase when ad-hoc measurements are not yet available.

Observing the small biases achieved by the microclimate system that embeds the PV model in Case 4 ($MBE = 0.5$ °C and $RMSE = 3.0$ °C), it is confirmed that the proposed ABL parameterization for the floating PV system is sufficiently efficient to assess the interaction of PV/atmosphere from the photovoltaic side point of view. Performance losses compared to Baseline 2 and Case 3 are associated with the computation of the local wind direction, friction velocity, and ambient temperature (given in the order of losses magnitude). Here, the single location for the large-scale meteorological measurement in B1 provides reliable data to represent the meteorology of the area, likely due to the representativeness of the site and the well-maintained sensors. However, determining the wind vector at the PV array level remains difficult: factors like sub-grid models (e.g., roughness lengths) and the quality of geometry and spatial discretization of detailed obstacles (e.g., hills, embankments) influence the calculation of the incoming wind velocity and airflow recirculation zones. In FPVs, the embankment's height relative to the water surface can change due to hydrological conditions or external activities like hydropower pumping, necessitating additional data to enhance the microclimate model's geometry and the local flow computation. Future improvements of the method should therefore focus on both the high-fidelity definition of local obstacles as well as the large-scale meteorological forcings by using spatially resolved sources like re-analysis dataset or multiple weather stations.
<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Evaluation of the model's capability in reproducing PV temperature for the module of reference. The boxplot indicates the distribution of module temperatures at the powerplant scale, computed by the microclimate model and using the proposed atmospheric parameterization.</figcaption></figure>

Case 4 also provides the spatial distribution of the PV module temperature $\partial_y T_p$ and Fig. 13 shows how the temperature distribution evolves throughout the day. The spreads in temperatures are found to be the largest during the peak hours of the sun, with a maximum magnitude of $5.8 \,^\circ\text{C km}^{-1}$ (i.e, inter-whisker interval). When both the irradiation level and the module temperature fall below $\Phi_{sw}^{\lambda P} < 500 \,\text{W m}^{-2}$ and $T_p < 50 \,^\circ\text{C}$, the distribution is narrowed to $1.3 \,^\circ\text{C km}^{-1}$ on average. Note that we focus solely on assessing the role of airflow with homogeneous radiative forcings at the level of the PV plant (see Fig. 9(a)). In fact during varying sky conditions, the effect of wind is likely to be little in front of the variation in the irradiation level; however, the computed temperature spread remains effective for clear-sky situations for which the usual big-PV assumption does not allow one to determine this information. With lowest module temperatures located around the array edge of atmospheric inflow, system-level modifications can be envisioned to reduce the effect from temperature mismatch on the electrical installation (e.g., inverter parameterization and PV string wiring). For the floating photovoltaic context, the spatial evolution in the PV temperature may also have a significant impact on the waterbody system, as the water temperature remains affected by the downward thermal radiation from the PV array, potentially promoting areas of reduced deficit in air-water temperatures.

### 3.4. Smooth-rough transitions reduce water savings from floating photovoltaics coverage

The evaporation of the waterbody covered by FPV panels is calculated using the relations Eqs. (3.4) and (3.6) using the microclimate model (case denoted 'C' hereafter), and we constructed a supplementary microclimate case without PV panels (uncovered case, U) to evaluate the impact of the coverage on total evaporation. Fig. 14 presents the evaporation level for both cases and for the day of interest. Surprisingly, only a slight decrease in total evaporation of 14% is due to solar panels, despite the water coverage reaching 66%. Note that this assessment excludes simulations conducted before 07:00 and after 16:00 ($\Phi_{sw}^{\lambda P} < 250 \,\text{W m}^{-2}$), avoiding timesteps where buoyancy-driven evaporation may be significant as the atmospheric parameterization for these conditions remains inadequate. Even if extending the analysis to the full dataset (e.g., including winter period with lower evaporation magnitudes) shows an increase from 14% to 20%, the results remain two times lower than the estimation of 42% provided in [46] for suspended coverage occupying 70% of the waterbody.

One reason for increased evaporation is that FPV absorb momentum energy and generate higher turbulence levels near the surface. As

depicted in Fig. 12, these levels vary but remain two to three times higher than those of free surfaces (i.e., Eq. (2.21)). The increase in shear amplifies vapour removal, especially in the surrounding areas of floating islands where no barriers for the transport of vapour limit the evaporation exchange. The effect of this smooth-rough transition is observed in Fig. 14 by looking at the evaporation level for the free water surfaces denoted 'FS' and the surface occupied by the PV panels 'PV', in both covered (C) and uncovered (U) cases. In particular, the free-water surface of the case with PV panels (C-FS) shows a surplus in evaporation of 8% to 10% on average compared to the free-water surface of the case without PV panels (U-FS). As the combination of higher shear stress with high mass-convection transfer rate depends on the PV array layout (e.g., island arrangement and sizes), it is likely that layouts with dispersed islands would obtain greater evaporation surplus. The FPV array remains a barrier for the transport of vapour from the waterbody to the atmosphere as seen by contrasting C-PV and U-PV, with a deficit in evaporation between 40% to 50% on average. In this case, the arrangement of the modules remains the principal factor driving the water saving potential (e.g. $H_p, \theta_p, S_p$). Finally, it should be highlighted that the lake temperature is imposed in the model, and further efforts should be made to implement a lake model to identify the influence of the water temperature on the evaporative performance. Knowing that the FPV system reduces the temperature of the water because of shadings, the 14% to 20% in total water savings may be greater in real-world applications.

## 4. Conclusion

This study presents a precursor CFD model for computing momentum absorption, convective heat transfer coefficients, and aerodynamic resistances in utility-scale PV arrays. Adapted for various tilt angles, module heights, and pitch distances, the method is particularly adapted to create wind-dependent CHTC functions for yield models and estimate the evaporation of partially covered waterbodies in floating PV context. We first applied the method to a wind-tunnel setup characteristic of a land-based application, demonstrating that altitude-based wind profile is accurately modelled above the array. Then reproducing the geometrical layout of a commercial FPV installation, we deduced the atmosphere boundary layer parameters adapted to shear dominated situations. Although a limited part of the numerically computed roughness length is cross-validated with a full-year set of turbulence measurements above an industrial FPV array, headwinds and tailwinds relative to the front surface of the PV modules are found to generate the
<figure><img src="image_18.png" /><figcaption>Fig. 14. Simulation results for the evaporation from covered and uncovered waterbodies, including total evaporation and sub-area evaporation per study cases.</figcaption></figure>

highest level of turbulence compared to sidewinds. The varying effect from wind-direction influences the local heat and moisture transfers by modulating the engine of the exchanges. Not yet studied in numerical weather modelling, it is also likely to change the atmosphere boundary layer in the vicinity of utility scale, promoting turbulent transfers in different regions around the PV array and based on the powerplant layout and the large-scale conditions dominating at the site.

The effect from wind directions is also explored using a microclimate CFD model that reproduces the local area of the industrial FPV array. Calculating the temperature of the PV panels and the evaporation of the waterbody, we estimated the large-scale module temperature gradient from 1.3 °C km⁻¹ on average to 5.8 °C km⁻¹, and a water saving rate between 14% and 20%, while the coverage is around two-thirds of the waterbody (clear-sky and sidewinds conditions). The upstream region is found to promote lower module temperatures whereas the downstream region remains at higher temperature levels. Measurements of module temperatures in the downstream region support this assertion, and further research is needed to observe this gradient on kilometre scales, as it could impact the performance of the photovoltaic asset (e.g., voltage imbalances). Finally, the spatial variation in smooth and rough surfaces, respectively free-water and PV surfaces) is found to increase the evaporation level of the free-water surfaces compared to situation of fully free-water surface by 8% to 10%. This effect, which does not cancel out the overall reduced evaporation due to the larger atmospheric resistance to vapour transport from the covered by PV areas, may become crucial if array layouts with disseminated floating islands are favoured.

### CRediT authorship contribution statement

**Baptiste Amiot:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. **Rémi Le Berre:** Writing – review & editing, Resources, Funding acquisition. **Stéphanie Giroux-Julien:** Supervision, Project administration, Conceptualization. **Martin Ferrend:** Writing – review & editing, Writing – original draft, Supervision, Software, Resources, Project administration, Methodology, Formal analysis, Conceptualization.

### Open access

For the purpose of Open Access, a CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

Writefull's model (including large-scale language-generation model and language checking) has been used to improve the writing style of this article. B. Amiot reviewed, edited, and revised the generated texts to his own liking and takes the ultimate responsibility for the content of this publication.

We thank G. Bayart and the EDF Power Solutions company for access to the operational data.

We thank B. Cuilhé (CEREA) for the support in running the microclimate digital twin.

### Annexes

#### A.1. Module temperature computation

In the microclimate model, a more precise description of the energy balance than Eq. (1.1) that includes energy lost in the infrared spectrum is used, it reads:

$$ \rho_p c_{p,p} \partial_t T_p = \alpha_{sw} \Phi_{sw}^{\text{p}} (1-\eta) - h_{cv,fr} (T_p - T_{c,l}) - h_{cv,re} (T_p - T_{c,l}) \\ - FF_{fr} \epsilon_{l,w} \sigma (T_p^4 - T_{sky}^4) - FF_{re} \epsilon_{l,w} \sigma (T_p^4 - T_{c,l}^4) \quad (A.1) $$

where $\rho_p$ is the PV module density, $c_{p,p}$ is the specific heat capacity of the PV module, $\alpha_{sw} = 0.9$ is the PV module absorptivity, $\Phi_{sw}^{\text{p}}$ is the downward solar radiation in the plane of array, $\eta = 0.185 (1 + 0.0028 \times (298.15 - T_p))$ is the module efficiency, $h_{cv,fr}$ and $h_{cv,re} = 1 \text{ W m}^{-2} \text{ K}^{-1}$ the convective heat transfer coefficients for the
front and the rear module faces respectively, $FF_{fr} = FF_{re} = 1$ are the form factors for the front and the rear face of the PV module, $T_{sky}$ is the sky temperature, $\epsilon_{lw} = 0.95$ is the PV module emissivity and $\sigma = 5.67 \times 10^{-8} W^2 m^{-1} K^{-4}$ is the Boltzmann constant. Assuming a steady state as well as a negligible relation between the module temperature and the electrical efficiency, the module temperature can be implicitly expressed so that it reads:

$$ (h_{cv,fr} + h_{cv,re}) T_p + (FF_{fr} + FF_{re}) \epsilon_{lw} \sigma T_p^4 = \alpha_{sw} \Phi_{sw}^{\lambda p} (1 - \eta) + (h_{cv,fr} + h_{cv,re}) T_{c,i} + FF_{fr} \epsilon_{lw} \sigma T_{sky}^4 + FF_{re} \epsilon_{lw} \sigma T_{c,i}^4 \quad (A.2) $$

For the purpose of the microclimate model, the radiation flux in the plane of array $\alpha_{sw} \Phi_{sw}^{\lambda p}$ and the downward infrared radiation $\Phi_{lw}^1 = \sigma T_{sky}^4$ are set constant and their values are those measured by the MS80S and SGR4 sensors, see Table 3, Figs. 8 and 9. Note that in Section 3.3, Baseline 1, Baseline 2 and Case 3 use Eq. (A.2) with $T_{c,i} = T_a$ (recorded at the meteostation level).

## A.2. Data processing for momentum

The wind turbulence varies significantly near the surface, making a single day's measurement, as shown in Fig. 9(b), inadequate to represent the roughness length of a given surface. To account for the photovoltaic array's effect, we utilized an expanded anemometer dataset covering March 2024 to March 2025. This resulted in more than 3.1 × 10⁸ observations for the wind components (bulk measurements) which are then postprocessed according to a set of quality criteria given as follows:

* $u_* > 0.1 \text{ m s}^{-1}$,

* $-100 < L_{MO} < 100$,

* $\sigma_{\theta_{lw}} < 45^{\circ}$,

where:

$$ u_* = \left( \overline{u_x' u_z'}^2 + \overline{u_y' u_z'}^2 \right)^{1/4} \quad (A.3) $$

given that $\overline{u_x' u_z'}$ and $\overline{u_y' u_z'}$ are the covariances of $u_x$, $u_z$ and $u_y$, $u_z$ that we measure over 10 min; $L_{MO}$ is the length of Monin-Obukhov computed from the bulk measurement of the wind components and the temperature deviations, and $\sigma_{\theta_{lw}}$ is the deviation from the mean wind direction during the period of integration. 7×10⁷ observations met these criteria which resulted in 12 389 datapoints averaged over 10 min periods. Roughness lengths are then discretized according to the wind direction. Each basket of wind direction covers a 30° range and is labelled with respect to the middle wind direction of the basket. In doing so, we obtain the roughness lengths distributions for the wind sectors that are deemed to be influenced by the PV array, they are shown in Fig. 11(a).

## Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## References

[1] D. Yang, W. Wang, C.A. Gueymard, T. Hong, J. Kleissl, J. Huang, M.J. Perez, R. Perez, J.M. Bright, X. Xia, D. van der Meer, I.M. Peters, A review of solar forecasting, its dependence on atmospheric sciences and implications for grid integration: Towards carbon neutrality, Renew. Sustain. Energy Rev. 161 (2022) 112348, http://dx.doi.org/10.1016/j.rser.2022.112348, URL https://www.sciencedirect.com/science/article/pii/S1364032122002593.

[2] C.F.M. Coimbra, Energy meteorology for the evaluation of solar farm thermal impacts on desert habitats, Adv. Atmospheric Sci. 42 (2) (2025) 313–326, http://dx.doi.org/10.1007/s00376-024-4242-3.

[3] C. Dupraz, H. Marrou, G. Talbot, L. Dufour, A. Nogier, Y. Ferard, Combining solar photovoltaic panels and food crops for optimising land use: Towards new agri-ovoltaic schemes, Renew. Energy 36 (10) (2011) 2725–2732, http://dx.doi.org/10.1016/j.renene.2011.03.005, URL https://www.sciencedirect.com/science/article/pii/S0960148111001194, Renewable Energy: Generation & Application.

[4] G. Tina, F. Bontempo Scavo, L. Merlo, F. Bizzarri, Analysis of water environment on the performances of floating photovoltaic plants, Renew. Energy 175 (2021) http://dx.doi.org/10.1016/j.renene.2021.04.082.

[5] A. Armstrong, N.J. Ostle, J. Whitaker, Solar park microclimate and vegetation management effects on grassland carbon cycling, Environ. Res. Lett. 11 (7) (2016) 074016, http://dx.doi.org/10.1088/1748-9326/11/7/074016.

[6] S. Yue, M. Guo, P. Zou, W. Wu, X. Zhou, Effects of photovoltaic panels on soil temperature and moisture in desert areas, Environ. Sci. Pollut. Res. 28 (14) (2021) 17506–17518, http://dx.doi.org/10.1007/s11356-020-11742-8, URL https://www.scopus.com/inward/record.uri?eid=2-s2.0-85099001000&doi=10.1007%2fs11356-020-11742-8&partnerID=40&md5=3f49a1822cf72223e13c4dc0bc4b3360, Cited by: 49.

[7] R.L. Nobre, C. Vagnon, S. Boulêtreau, F. Colas, F. Azémar, L. Tudesque, N. Parthuisot, P. Millet, J. Cucherousset, Floating photovoltaics strongly reduce water temperature: A whole-lake experiment, J. Environ. Manag. 375 (2025) 124230, http://dx.doi.org/10.1016/j.jenvman.2025.124230, URL https://www.sciencedirect.com/science/article/pii/S0301479725002063.

[8] G. Exley, R. Hernandez, T. Page, M. Chipps, S. Gambro, M. Hersey, R. Lake, K.-S. Zoannou, A. Armstrong, Scientific and stakeholder evidence-based assessment: Ecosystem response to floating solar photovoltaics and implications for sustainability, Renew. Sustain. Energy Rev. 152 (2021) 111639, http://dx.doi.org/10.1016/j.rser.2021.111639, URL https://www.sciencedirect.com/science/article/pii/S136403212100914X.

[9] G.A. Barron-Gafford, R.L. Minor, N.A. Allen, A.D. Cronin, A.E. Brooks, M.A. Pavao-Zuckerman, The photovoltaic heat island effect: Larger solar power plants increase local temperatures, Sci. Rep. 6 (2016) http://dx.doi.org/10.1038/srep35070.

[10] A.M. Broadbent, E.S. Krayenhoff, M. Georgescu, D.J. Sailor, The observed effects of utility-scale photovoltaics on near-surface air temperature and energy balance, J. Appl. Meteorol. Climatol. 58 (5) (2019) 989–1006, http://dx.doi.org/10.1175/jamc-d-18-0271.1.

[11] P. Li, X. Gao, Z. Li, T. Ye, X. Zhou, Effects of fishery complementary photovoltaic power plant on radiation, energy flux and driving forces under different synoptic conditions, Sci. Rep. 13 (1) (2023) 9084, http://dx.doi.org/10.1038/s41598-023-36314-x.

[12] I. Cesar, B. Van Aken, Evaluation method and module design for cost-effective compliance with irradiance guidelines to maintain soil quality in solar parks, EPJ Photovolt. 16 (2025) http://dx.doi.org/10.1051/epjpv/2025003.

[13] L. Guoqing, R.R. Hernandez, G.A. Blackburn, G. Davies, M. Hunt, J.D. Whyatt, A. Armstrong, Ground-mounted photovoltaic solar parks promote land surface cool islands in arid ecosystems, Renew. Sustain. Energy Transition 1 (2021) http://dx.doi.org/10.1016/j.rset.2021.100008, URL https://www.scopus.com/inward/record.uri?eid=2-s2.0-85138594845&doi=10.1016%2fj.rset.2021.100008&partnerID=40&md5=78a7c67bd025886db9be348551291dce, Cited by: 13; All Open Access, Gold Open Access.

[14] Y. Hu, C. Su, Y. Zhang, Research on the microclimate characteristics and cold island effect over a reservoir in the hexi region, Adv. Atmospheric Sci. 5 (1) (1988) 117–126, http://dx.doi.org/10.1007/BF02657351, URL https://www.scopus.com/inward/record.uri?eid=2-s2.0-0024253322&doi=10.1007%2FBF02657351&partnerID=40&md5=154dab1f87694b2875d6e9de3126b2bc, Cited by: 23.

[15] A. Driesse, M. Theristis, J.S. Stein, PV module operating temperature model equivalence and parameter translation, 2022, http://dx.doi.org/10.1109/PVSC48317.2022.9938895.

[16] M.G. Farr, J.S. Stein, Spatial variations in temperature across a photovoltaic array, 2014, http://dx.doi.org/10.1109/PVSC.2014.6925301.

[17] M.M. Escribano, M.G. Solano, Í. de la Parra Laita, J.M. Álvarez, L. Marroyo, E.L. Piguerales, Module temperature dispersion within a large PV array: Observations at the Amareleja PV plant, IEEE J. Photovolt. 8 (2018) 1725–1731, URL https://apisemanticscholar.org/CorpusID:53095496.

[18] M. Maeda, A. Nagaoka, K. Araki, K. Nishioka, Practical and simplified measurements for representative photovoltaic array temperatures robust to climate variations, Sol. Energy 231 (October 2020) (2022) 243–251, http://dx.doi.org/10.1016/j.solener.2021.11.035.

[19] H. Goverde, D. Goossens, J. Govaerts, F. Catthoor, K. Baert, J. Poortmans, J. Driesen, Spatial and temporal analysis of wind effects on PV modules: Consequences for electrical power evaluation, Sol. Energy 147 (2017) 292–299, http://dx.doi.org/10.1016/j.solener.2016.12.002, URL https://www.scopus.com/inward/record.uri?eid=2-s2.0-850160049188&doi=10.1016%2fj.solener.2016.12.002&partnerID=40&md5=db09a546aef275d04499115493440279, Cited by: 34.

[20] F.M. Daniel-Durandt, A.J. Rix, The long-term influence of wind and temperature on performance and degradation within an utility-scale photovoltaic plant, IET Renew. Power Gener. 19 (e13182) (2025) http://dx.doi.org/10.1049/rpg2.13182, Published online: 2025.
[21] E. Skoplaki, J.A. Palyvos, Operating temperature of photovoltaic modules: A survey of pertinent correlations, Renew. Energy 34 (1) (2009) 23–29, http://dx.doi.org/10.1016/j.renene.2008.04.009.

[22] V. Masson, M. Bonhomme, J.-L. Salagnac, X. Briottet, A. Lemonsu, Solar panels reduce both global warming and urban heat island, Front. Environ. Sci. 2 (JUN) (2014) http://dx.doi.org/10.3389/fenvs.2014.00014.

[23] H. Taha, Episodic performance and sensitivity of the urbanized MM5 (uMM5) to perturbations in surface properties in houston texas, Bound.-Layer Meteorol. 127 (2) (2008) 193–218, http://dx.doi.org/10.1007/s10546-007-9258-6.

[24] R. Macdonald, R. Griffiths, D. Hall, An improved method for the estimation of surface roughness of obstacle arrays, Atmos. Environ. 32 (11) (1998) 1857–1864, http://dx.doi.org/10.1016/S1352-2310(97)00403-2, URL https://www.sciencedirect.com/science/article/pii/S1352231097004032.

[25] R. Chang, Y. Yan, Y. Luo, C. Xiao, C. Wu, J. Jiang, W. Shi, A coupled WRF-PV mesoscale model simulating the near-surface climate of utility-scale photovoltaic plants, Sol. Energy 245 (2022) 278–289, http://dx.doi.org/10.1016/j.solener.2022.09.023.

[26] A.M. Aly, On the evaluation of wind loads on solar panels: The scale issue, Sol. Energy 135 (2016) 423–434, http://dx.doi.org/10.1016/j.solener.2016.06.018, URL https://www.sciencedirect.com/science/article/pii/S0038092X16301979.

[27] C.M. Jubayer, K. Siddiqui, H. Hangan, CFD analysis of convective heat transfer from ground mounted solar panels, Sol. Energy 133 (2016) 556–566, http://dx.doi.org/10.1016/j.solener.2016.04.027.

[28] V. Fthenakis, Y. Yu, Analysis of the potential for a heat island effect in large solar farms, 1 (1) (2013) 3362–3366. http://dx.doi.org/10.1109/PVSC.2013.6745171.

[29] D. Lindholm, J. Selj, T. Kjeldstad, H. Fjær, V. Nysted, CFD modelling to derive U-values for floating PV technologies with large water footprint, Sol. Energy 238 (April) (2022) 238–247, http://dx.doi.org/10.1016/j.solener.2022.04.028.

[30] B.J. Stanislawski, T. Harman, T.J. Silverman, R.B. Cal, M. Calaf, Row spacing as a controller of solar module temperature and power output in solar farms, J. Renew. Sustain. Energy 14 (6) (2022) 063702, http://dx.doi.org/10.1063/5.0123641, arXiv:https://pubs.aip.org/aip/jrse/article-pdf/doi/10.1063/5.0123641/16592541/063702_1_online.pdf.

[31] S.E. Smith, B.J. Stanislawski, B.K. Eng, N. Ali, T.J. Silverman, M. Calaf, R.B. Cal, Viewing convection as a solar farm phenomenon broadens modern power predictions for solar photovoltaics, J. Renew. Sustain. Energy (2022) URL https://apisemanticscholar.org/CorpusID:254357312.

[32] A. Glick, N. Ali, J. Bossuyt, G. Recktenwald, M. Calaf, R.B. Cal, Infinite photovoltaic solar arrays: Considering flux of momentum and heat transfer, Renew. Energy 156 (2020) 791–803, http://dx.doi.org/10.1016/j.renene.2020.03.183.

[33] A. Glick, S.E. Smith, N. Ali, J. Bossuyt, G. Recktenwald, M. Calaf, R.B. Cal, Influence of flow direction and turbulence intensity on heat transfer of utility-scale photovoltaic solar farms, Sol. Energy 207 (November 2019) (2020) 173–182, http://dx.doi.org/10.1016/j.solener.2020.05.061.

[34] code_saturne official website, 2025, URL https://www.code-saturne.org/cms/web/.

[35] V. Guimet, D. Laurence, A linearised turbulent production in the $k-\epsilon$ model for engineering applications, in: Engineering Turbulence Modelling and Experiments 5, Elsevier, 2002, pp. 157–166.

[36] F.R. Menter, Two-equation eddy-viscosity turbulence models for engineering applications, AIAA J. 32 (1994) 1598–1605, URL https://api.semanticscholar.org/CorpusID:120712103.

[37] H. Amino, C. Flageul, S. Benhamadouche, I. Tiselj, B. Carissimo, M. Ferrand, A time-staggered second order conservative time scheme for variable density flow, Internat. J. Numer. Methods Fluids 94 (12) (2022) 1964–1995, http://dx.doi.org/10.1002/fld.5116.

[38] Salome official website, 2025, URL https://www.salome-platform.org/?lang=en.

[39] J. Silva, C. Ribeiro, R. Guedes, M.-C. Rua, F. Ulrich, Roughness length classification of Corine Land Cover classes, in: Proceedings of EWEC 2007, 2007.

[40] A.K. Betts, Non-precipitating cumulus convection and its parameterization, Q. J. R. Meteorol. Soc. 99 (419) (1973) 178–196, http://dx.doi.org/10.1002/qj.49709941915.

[41] E. Bouzereau, Représentation Des Nuages Chauds Dans Le Modèle Météorologique “Mercure” : Application Aux Panaches D’aéroréfrigérants Et Aux Précipitations Orographiques (Theses), Université Pierre et Marie Curie - Paris VI, 2004, URL https://theses.hal.science/tel-00009305.

[42] H. Charnock, Wind stress on a water surface, Q. J. R. Meteorol. Soc. 81 (350) (1955) 639–640, http://dx.doi.org/10.1002/qj.49708135027.

[43] D.J.G.T.O. Glanz, Lincoln Lake ecologic study, Water Resour. Eng. (1973).

[44] F.L. Test, R.C. Lessmann, A. Johary, Heat transfer during wind flow over rectangular bodies in the natural environment, J. Heat Transfer 103 (2) (1981) 262–267, http://dx.doi.org/10.1115/1.3244451.

[45] I. Peters, A. Nobre, Deciphering the thermal behavior of floating photovoltaic installations, Sol. Energy Adv. 2 (November 2021) (2022) 100007, http://dx.doi.org/10.1016/j.seja.2021.100007.

[46] F. Bontempo Scavo, G.M. Tina, A. Gagliano, S. Nižetić, An assessment study of evaporation rate models on a water basin with floating photovoltaic plants, Int. J. Energy Res. 45 (1) (2021) 167–188, http://dx.doi.org/10.1002/er.5170, URL https://onlinelibrary.wiley.com/doi/10.1002/er.5170.