Assimilation of temperature and relative humidity observations from personal weather stations in AROME-France

Personal weather station (PWS) networks owned by citizens now provide near-surface observations at a spatial density unattainable with standard weather stations (SWSs) deployed by national meteorological services. This article aims to assess the benefits of assimilating PWS observations of screen-level temperature and relative humidity in the AROME-France model in the same framework of experiments carried out to assimilate PWS observations of surface pressure in a previous work. Several methods for pre-processing these observations, in addition to the usual data assimilation (DA) screening, are evaluated and selected. After pre-processing, nearly 4700 PWS temperature and 4200 PWS relative humidity observations are assimilated per hour, representing 3 and 6 times more than SWS observations, respectively. Separate assimilation of each variable in the atmosphere with the three-dimensional ensemble variational (3DEnVar) DA scheme significantly reduces the root-mean-square deviation between SWS observations and forecasts of the assimilated variable at 2 m height above ground level up to 3 h of forecasts. Improvements to the near-surface temperature and relative humidity fields analysed are shown for a sea breeze case during a heatwave and a fog episode. However, degradation of short-range forecasts are found when PWS observations are assimilated with the current operational 3DVar DA scheme in the atmosphere or jointly in the atmosphere and at the surface with 3DEnVar and optimal interpolation DA schemes. These results demonstrate that the benefit of assimilating PWS temperature and relative humidity observations can be highly dependent on the DA schemes and settings employed.

1 Introduction

The increase in spatial and temporal resolution of regional numerical weather prediction (NWP) models requires their analyses to be initialized by spatially and temporally dense observations to represent meteorological phenomena on increasingly finer scales (Gustafsson et al., 2018). At Météo-France, the operational regional NWP system designed to forecast such phenomena with an up to 51 h lead time is AROME-France (Seity et al., 2011; Brousseau et al., 2016). AROME-France currently uses a three-dimensional variational (3DVar) DA scheme. This scheme is planned to be replaced by a three-dimensional ensemble variational (3DEnVar) DA scheme, with an up-to-date background error covariance matrix, improving the spread of the information from the observations (Montmerle et al., 2018; Michel and Brousseau, 2021).

Phenomena whose representation in the analyses needs to be improved are meso-γ- to meso-β-scale phenomena (2 to 200 km; Orlanski, 1975) such as thunderstorms, breezes or fog (Stull, 1988) which cause substantial thermodynamic changes, particularly in the atmospheric boundary layer. New observations are emerging to improve analyses above the surface, from ground-based microwave radiometers (Caumont et al., 2016; Bell et al., 2022; Vural et al., 2024), water vapour lidars (Flamant et al., 2021), or aircraft (Pourret et al., 2022). Even if these instruments provide dense vertical observations, their horizontal density remains low. Work is underway to exploit horizontally dense satellite observations near the surface, such as radiances from the Meteosat Spinning Enhanced Visible and Infrared Imager (SEVIRI) (Sassi et al., 2019), which are currently assimilated with a spacing of the order of 30 km and are still difficult to assimilate at full resolution.

The growth of crowdsourced data, i.e. the data obtained from a range of sensors belonging to the public and shared via the Internet, could play a significant role in improving near-surface analyses of NWP models, in particular data from personal weather stations (PWSs) (Hintz et al., 2019; Coney et al., 2022). In order to obtain observations comparable with standard weather station (SWS) observations, PWS observations should be properly pre-processed, i.e. bias-corrected or quality-controlled or both, the bias designating the deviation of a PWS observation time series from a reference time series (other observations or a modelled field). Indeed, PWSs provide observations that may not comply with World Meteorological Organization (WMO) methods: PWS sensors have heterogeneous environments of siting, have sometimes unsuitable shelters, can be positioned at various heights above ground level (a.g.l.), and can suffer from time lag compared to sensors following WMO requirements (Bell et al., 2015; Varentsov et al., 2020; Fenner et al., 2021). After pre-processing, observational studies showed that PWS observations, when combined with SWS observations, are able to describe near-surface thermodynamic variations associated with mesoscale phenomena such as thunderstorms that are partially visible with SWS observations only (Clark et al., 2018; Mandement and Caumont, 2020).

Various bias correction (BC) methods have been developed. The bias can be computed over a short period (e.g. 6 h for Clark et al., 2018, or approximately 24 h for Mandement and Caumont, 2020). To account for the effects of an inappropriate sheltering or siting, the bias can be computed as a linear function of solar incident radiation or by a more sophisticated decomposition, e.g. with multilinear functions (Sgoff et al., 2022), generalized additive mixed models (Cornes et al., 2019), or machine learning methods (Beele et al., 2022; Marquès, 2023).

Regarding quality control (QC) methods, CrowdQC (Meier et al., 2017) and CrowdQC+ (Fenner et al., 2021) were designed to remove PWS temperature observations deemed erroneous based on their deviation from neighbouring observations. CrowdQC was used in numerous studies of urban temperatures (Feichtinger et al., 2020; Venter et al., 2020; Potgieter et al., 2021; Zumwald et al., 2021). Another QC for PWS temperature observations, used operationally in post-processing algorithms by MET Norway, is called Titan (QC-Titan hereafter; Båserud et al., 2020; Nipen et al., 2020). CrowdQC+ and QC-Titan are comparable in a way as they both use the information from their neighbours to remove inconsistent observations, e.g. the m5 procedure of CrowdQC+ is comparable to the spatial buddy check of QC-Titan (Fenner et al., 2021). In addition, as CrowdQC+ uses local climate zones, it is challenging to use on a scale larger than a city. For both PWS temperature and relative humidity observations, Mandement and Caumont (2020) have developed a QC (QC-MC hereafter) using an adaptive rejection threshold based on a comparison to interpolated SWS observations. Both QC-Titan and QC-MC are tested in this study.

Once PWS observations have been pre-processed, the opportunity of their assimilation arises. Sgoff et al. (2022) carried out experiments assimilating pre-processed (bias-corrected and quality-controlled) PWS temperature and relative humidity observations, using the flow-dependent local ensemble transform Kalman filter DA scheme of the Icosahedral Nonhydrostatic model with 2 km resolution (ICON-D2). They showed that the assimilation of PWS observations is beneficial; however, no experiments assimilating at the same time both SWS and PWS observations were conducted. On the other hand, Demortier et al. (2024, hereafter D24) carried out experiments assimilating simultaneously pre-processed PWS and SWS observations of surface pressure with both 3DVar and 3DEnVar DA schemes in AROME-France. Statistically significant improvements in mean sea level pressure forecasts of up to 9 h range were found with the 3DEnVar scheme.

This article extends the Sgoff et al. (2022) and D24 works by simultaneously assimilating pre-processed PWS and SWS observations of screen-level temperature and relative humidity and addresses the following questions. What are the most effective pre-processing methods for assimilating PWS data? Are AROME-France's analyses and forecasts improved by assimilating pre-processed PWS temperature and relative humidity observations? What impact does the choice of the DA scheme have?

The remainder of this article is organized as follows. Section 2 describes the AROME-France NWP system and its DA schemes in the atmosphere and at the surface, as well as the observations used and the pre-processing methods of these observations. Section 3 describes the assimilation experiments. Objective results of these experiments are given in Sect. 4, and results for case studies are detailed in Sect. 5. Finally, results are summarized and discussed in Sect. 6.

2 Datasets and methods

The assimilation experiments run from 6 September to 5 October 2021, as in D24. This 1-month study period encompasses a diverse range of meteorological events. It includes two periods of anticyclonic conditions, during which a heatwave (from 2 to 8 September) and fog (from 22 to 24 September) were observed. It also includes disturbed weather episodes such as squall lines on 8 September in the south-west of France, a mesoscale convective system in the south-east of France on 14 September and heavy precipitation events from 2 to 4 October.

2.1 AROME-France NWP system

AROME-France is the limited-area NWP model developed by Météo-France, operational since December 2008 (Seity et al., 2011; Brousseau et al., 2016). AROME-France is a spectral model coupled to the global NWP system ARPEGE (Courtier et al., 1998). AROME-France (AROME hereafter) has a 1.3 km grid spacing on the horizontal and has 90 vertical levels in the atmosphere ranging from 5 m height a.g.l. (lowest atmospheric model level) up to 10 hPa.

AROME uses a mass-flux shallow convection scheme (Pergaud et al., 2009). Its associated surface scheme is SURFEX (Masson et al., 2013). AROME and SURFEX exchange flux and near-surface variables such as temperature, humidity, and wind. Surface variables are computed over tiles representing four types of surfaces: nature, town, inland water, and ocean. The surface layers have a dedicated assimilation system, which is described in Sect. 2.1.2.

AROME has a 1 h DA cycle for the atmosphere and a 3 h DA cycle for the surface. From the analysis, the experiments in this study produce 24 h forecasts at 00:00, 06:00, 12:00, and 18:00 UTC.

In AROME, an average of 49 253 observations (active; see below) are assimilated per hour during the 1-month study period. A total of five observation systems account for 99 % of incoming active observations: radar (56.5 %), SWSs (17.0 %), satellites (9.8 %), radiosoundings (9.6 %), and aircraft (6.1 %). Among the variables observed by SWSs, the ones assimilated in AROME are screen-level temperature and relative humidity, zonal and meridional wind at 10 m height a.g.l., and the surface geopotential. More details are given in Sect. 2.2.1.

2.1.1 AROME atmospheric DA system

The atmospheric DA procedure is composed of two main steps, which are the screening, which is a series of QC checks, and the minimization of the cost function used for variational DA (see D24). During both steps, modelled variables are compared to observed variables by an observation operator. AROME has 13 prognostic variables and 5 control variables, which include temperature and specific humidity. Relative humidity is a diagnostic variable computed from specific humidity, temperature, and pressure. The formula used in AROME is given in Appendix A. The model equivalents of observed screen-level temperature (T) and relative humidity (RH) are 2 m height a.g.l. temperature (T_2 m) and relative humidity (RH_2 m), described by Vasiljevic et al. (1992) and ECMWF (2023). The computation of T_2 m uses the temperature from the lowest atmospheric model level and the surface temperature, based on Monin–Obukhov theory – with separate solutions for stable and unstable conditions (Cardinali et al., 1994). RH_2 m is computed using only the lowest level fields of the atmospheric model. Thus, the effect of the type of surface on observations is taken into account by these observation operators, either directly when it includes a surface variable modelled by SURFEX (for T_2 m) or indirectly when it only includes atmospheric variables, modified through the surface–atmosphere exchanges during the forecast (for RH_2 m).

The difference between the observation and the model equivalent from the background (1 h forecast starting from the previous analysis), i.e. the observation minus the background, is further called OmB. In the same way, the difference between the observation and its corresponding value in the analysis is further called OmA.

One of the main QC checks of the screening rejects observations, verifying

$$\begin{array}{}\text{(1)}& \Vert \text{OmB}\Vert >{\displaystyle \frac{\mathit{\alpha }}{{\mathit{\sigma }}_{\text{coef}}}}\sqrt{{\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{2}}},\end{array}$$

where α=4 and σ_coef=0.9 are two coefficients, σ_o is the standard deviation of observation errors set to 1.4 °C and 10 % for screen-level T and RH observations, σ_b is the standard deviation of background errors, and the right-hand term is referred to as the rejection threshold. As an example, for T_2 m, during the first hour of the assimilation experiments (i.e. 6 September 2021 at 01:00 UTC), the average σ_b equals 0.5 °C and the rejection threshold is up to 6.5 °C. For RH_2 m, the average σ_b equals 10 % and the rejection threshold is up to 70 %. These rejection thresholds are high because current near-surface assimilated observations (SWS observations) are considered as anchor observations for the model, and screening has been designed to keep as much of these observations as possible.

All active observations, i.e. those that have passed the screening, are combined with the background to produce the analysed state of the atmosphere. AROME currently uses a 3DVar DA scheme to constitute the atmospheric analysis fields. The only difference between the 3DVar and the 3DEnVar DA schemes comes from the background error covariance matrix B, which is static (prescribed) for the 3DVar and dynamic for the 3DEnVar, estimated at each hour from the AROME Ensemble Data Assimilation system. To illustrate the difference between the two schemes, the increments (i.e. the analysis minus the background) at the lowest atmospheric model level when single screen-level T and RH observations are assimilated in two idealized experiments are shown in Fig. 1. The propagation of the information from an observation is isotropic with the 3DVar DA scheme (Fig. 1a and c), whereas with the 3DEnVar DA scheme it is anisotropic and localized within 25 km around the observation (Fig. 1b and d). At the observation location, the values of the increments differ. From an identical OmB of 3.4 °C in both experiments, the increment at the observation location reaches 1 °C with the 3DVar DA scheme and is reduced to 0.4 °C with the 3DEnVar DA scheme. Similarly, the assimilation of a single RH_2 m observation results in an increment of 3.3 % (1.1 %) with the 3DVar (3DEnVar) DA scheme. The relative observation weight is diminished from the 3DVar to the 3DEnVar DA scheme due to lower σ_b on average during the study period for temperature observations. Further details on the 3DVar and the 3DEnVar DA scheme configurations are provided by D24.

Figure 1 Increments at the lowest atmospheric model level of (a, b) temperature and (c, d) specific humidity, resulting from the assimilation of screen-level T and RH observations, respectively, at Blagnac SWS on 6 September 2021 at 01:00 UTC with the (a, c) 3DVar and (b, d) 3DEnVar DA schemes. The black circle indicates a 25 km distance around the observation location.

2.1.2 AROME surface DA system

AROME has a dedicated surface DA scheme for the initialization of the soil prognostic variables, which consist of two ground temperature variables and two groundwater content variables (Giard and Bazile, 2000). The surface DA scheme has its own screening, which uses the same formulation as the atmospheric one (Eq. 1), with, for T_2 m (RH_2 m), α set to 5 (2.5), σ_coef set to 1, σ_o set to 1.3 °C (10 %), and σ_b set to 1.6 °C (10 %).

The surface DA scheme is based on atmospheric increments of T_2 m and RH_2 m. First, T_2 m and RH_2 m are computed with the observation operator described in Sect. 2.1.1. Then, a univariate optimal interpolation (OI) scheme combines screen-level observations with T_2 m and RH_2 m background, using isotropic structure functions. The characteristic distance D used in the OI scheme is set to 100 km for T_2 m and RH_2 m observations, which is a compromise between the resolution of the NWP model, the scale of the weather phenomenon represented and the density of the weather stations. Then T_2 m and RH_2 m increments are used to correct the soil variables with linear interpolations. More details on the surface DA are given by Sassi et al. (2019).

2.1.3 A posteriori diagnostics

A number of diagnostic tools exist to understand a posteriori the way in which the DA schemes use the observations. Both the Desroziers diagnostic of observation error and the spatial Desroziers diagnostic (also used in D24) are used (Desroziers et al., 2005). The spatial Desroziers diagnostic is computed to define a minimum thinning length λ from which the correlated observation errors are found to be significantly low. Then, the observation network is thinned by selecting one random observation per mesh from a horizontally spaced λ grid (D24).

2.2 Screen-level observations of temperature and relative humidity

A distinction is made here between standard weather stations and personal weather stations: the former are bought, operated and maintained by national meteorological and hydrological services, while the latter are managed by private individuals or third-party organizations. Behind this distinction, it is important to keep in mind that citizens or third-party organizations can acquire stations identical to meteorological services, just as meteorological services can use low-cost stations dedicated to citizens.

2.2.1 Standard weather station observations

This study uses only SWS observations assimilated in AROME, which constitute a subset of all SWS observations. On average over 1 month, the number of SWSs providing at least one observation per hour of temperature and relative humidity is 2440 and 1600, respectively. Only 65 % of SWS temperature observations entering the AROME DA system are located in France. This is only 48 % for relative humidity. Over France, there are more than twice as many observations of temperature as there are observations of relative humidity, and both are fairly evenly distributed (Fig. 2). Regarding the height of observation, air temperature and relative humidity should be and are generally measured between 1.25 and 2 m height a.g.l. (World Meteorological Organization, 2023), the rule being 1.5 m in France, but often a little higher up in mountainous areas with heavy snowfall. For SWS wind speed observations used in Sect. 4, wind speed is generally measured at 10 m height a.g.l.

Figure 2 Number of (a) SWSs measuring temperature, (b) SWSs measuring relative humidity, (c) Netatmo PWSs, and (d) StatIC PWSs providing at least one observation during the 1-month study period over France. Observation counts are binned into approximately 0.15 × 0.1° bins. The dashed rectangles delimit the four domains used in the cases studied: (1) Paris, (2) Toulouse, (3) along the Atlantic coast for the heatwave case, and (4) the Saône Valley for the fog case.

2.2.2 Personal weather station observations

Three networks of PWSs are used in this study: Netatmo PWSs, StatIC PWSs, and Toulouse Métropole PWSs. In contrast to SWS observations, PWS observations are not currently assimilated operationally in AROME.

The Netatmo PWS network is the largest available in near-real time in France (Mandement and Caumont, 2020). As they are owned by citizens, PWSs are unevenly distributed in France (Fig. 2). The outdoor module of the Netatmo PWS network contains temperature and relative humidity sensors. On average over 1 month, the number of PWSs providing at least one observation per hour of temperature and relative humidity is 63 099 and 63 117, respectively. With appropriate sheltering, these sensors have a median and a 95 % range of departures from a reference sensor of about 0 °C ± 0.9 °C in temperature and 3 % ± 7 % in relative humidity, showing their correct intrinsic quality (Mandement and Caumont, 2020). However, the stations as they are sold are not sheltered according to WMO recommendations: short-wave and long-wave radiation affect T and RH, causing departures from sheltered reference sensors (Büchau, 2018). For temperature, these departures are generally overestimations of warm temperatures and underestimations of cold temperatures. These departures are complex depending on the location of the station in relation to its close environment (e.g. wall, balcony, garden; Varentsov et al., 2020) and its actual height above ground level.

The StatIC PWS network from the Infoclimat association encompassed 880 stations in 2023 (Garcelon et al., 2023), almost 600 of which are free to download over the period via their open-data portal. It gathers Davis Instruments Vantage Pro2 and Vantage Vue PWSs, as well as sensors from the companies Dragino or Talkpool in shelters, owned by citizens or by the association itself. To join the network, each station should follow the WMO recommendations adapted by Météo-France (Leroy, 2014) wherever possible: T and RH observations in a normalized shelter, in a clear environment, and at 1.5 m height a.g.l. However, in enclosed urban environments (with low air circulation), the recommended height is 1.5 m above the roof, i.e. various heights a.g.l. depending on the height of the roof. StatIC PWSs are set up in meteorologically interesting locations such as La Chapelle-en-Vercors, which is a cold hole, or areas not well covered by SWSs. The accuracy given for the Davis Instruments Vantage Pro2 or Vue PWS is 0.3 °C for temperature and 2 % for relative humidity (Garcelon et al., 2023).

Toulouse Métropole has deployed Davis Vantage Pro2 PWSs around the city of Toulouse for the study of the urban heat island (Dumas et al., 2021). During the study period, 33 PWSs are available.

Figure 3 Histograms of average OmB over the 1-month study period for (a, d) Netatmo PWS, (b, e) StatIC PWS, and (c, f) SWS observations of (a–c) temperature and (d–f) relative humidity. Boxes indicate the number of OmB time series, μ the average OmB for all time series, and σ the standard deviation of average OmB for all time series.

Figure 3 shows the average OmB distributions for the two largest PWS networks and the SWS network. For temperature over the study period, there are only one-fifth as many StatIC PWS observations as SWS observations, while there are about 25 times as many Netatmo PWS observations. However, the heterogeneous siting of the PWSs, seen by the OmB standard deviation, increases from SWSs to StatIC PWSs to Netatmo PWSs. The results are similar for relative humidity.

In this study, the only PWS observations assimilated are from Netatmo PWSs, hereafter referred to as PWSs. StatIC and Toulouse Métropole PWS observations are used to evaluate pre-processing methods of PWS observations of temperature. For the evaluation, Netatmo PWS observations are interpolated at the location of StatIC and Toulouse Métropole PWSs, providing an estimate which is compared to the StatIC and Toulouse Métropole PWS observations. The linear interpolation method takes into account the vertical profile of temperature. The use of independent (i.e. non-assimilated) PWS observations allows the absence of bias in the model to be verified.

2.3 Satellite observations

Observations of top-of-the-atmosphere bidirectional reflectance from the high-resolution visible (HRV) channel from the SEVIRI radiometer, on board the Meteosat Second Generation 3 (MSG-3) satellite (also called Meteosat-10) in rapid scanning service at 1 km horizontal resolution every 5 min, are used in Sect. 5.2 to locate clouds. The HRV channel is a broadband channel which is sensitive to wavelengths between 0.4 and 1.1 µm (Schmetz et al., 2002). Reflectance observations from EUMETSAT are corrected depending on the solar angle following the method of Li and Shibata (2006) applied to visible channels in Météo-France operational products (Derrien and Le Gléau, 2010).

An estimate of the global solar radiation at each PWS location has been derived from the spatialized global solar radiation product made by Météo-France in near-real time. This product uses both in situ observations and surface solar irradiance from MSG satellites (EUMETSAT, 2017). Nighttime values are set to zero. The hourly global solar irradiance is cumulated over each hour and is used as a predictor for one of the bias correction methods of PWS observations.

2.4 Pre-processing methods

In the article, PWS observations used are interpolated to round hours, and PWSs with identical coordinates are removed, as in the step preparation of D24. Hereafter, these observations, prepared but not pre-processed, are referred to as raw PWS observations.

2.4.1 Bias correction (BC) methods

The bias correction (BC) ensures that possible biases (i.e. a systematic departure from a reference) in PWS observation time series do not propagate to the model. In this study, raw PWS observations of screen-level temperature and relative humidity have been monitored – i.e. T_2 m and RH_2 m OmB time series are computed – in the operational AROME 3DVar configuration, which is called the Monitor experiment (as in Sgoff et al., 2022, and D24). Four types of BC methods that can be applied in real time because they use past OmB statistics, are tested. They are named as follows:

• BC-D removes the average OmB of the last 24 h (average of 24-hourly OmBs, D for daily), in a manner close to Mandement and Caumont (2020).

• BC-M removes the average OmB of the last 30 d (M for monthly).

• BC-S uses the BC method developed by Sgoff et al. (2022), in which S signifies Sgoff. The bias is computed from five predictors composed of trigonometric basis functions from the hour of the measured observations with 2 d inertia.

• BC-R uses a random forest method trained with the hour of the day and cumulative hourly global solar irradiance predictors from the last 30 d, close to the BC used by Beele et al. (2022) but with fewer predictors (R for radiation).

2.4.2 Quality control (QC) methods

Quality control (QC) methods are used in this study to remove PWS observations judged erroneous because the current screening for near-surface observations has not been designed to deal with observations with OmB statistics so different from those of SWS (Sect. 2.2.2).

Titan or titanlib (Båserud et al., 2020) is an open-source library composed of various checks or tests. In this article, as in Nipen et al. (2020), Titan-QC designates the combination of three main spatial tests, which are the buddy check, the spatial consistency test, and the isolation test. The buddy check consists of removing observations if their deviation from the average is more than twice the standard deviation of the observations in the neighbourhood within a 15 km radius. The spatial consistency test consists of an iterative cross-validation procedure. For each observation, it adjusts a vertical profile, and computes an estimation of a true observation, and the standard deviation of the observations in the neighbourhood. The observation is removed if the ratio of the squared deviation from the observation estimation divided by the standard deviation is greater than 4 (or 8 for negative values). The process is repeated until no observations are removed. The isolation test removes observations which have fewer than five stations within a 15 km radius and 200 m elevation difference.

QC-MC (Mandement and Caumont, 2020) has been adapted for use on a larger scale, throughout France even though it was designed for the scale of a French region: computations are done simultaneously for 11 selected climatological areas. The second main adaptation is the use, for each PWS observation, of OmB statistics from AROME instead of comparisons between this PWS observation and an estimate derived from neighbouring SWS observations. In QC-MC, every hour, in every climatological area, for each PWS, the rms's (root mean squares) of PWS OmBs are computed on the last 6 h to separate observations that differ from the background at a given time due to physical reasons (e.g. a phenomenon not well positioned in the background) from the observations that differ all the time from the background. PWSs having rms OmB exceeding an adaptive threshold, computed according to Mandement and Caumont (2020), are removed.

3 Assimilation experiment design

3.1 Choice of bias correction and quality control

In the literature, some authors carry out a QC without a BC (Nipen et al., 2020). To evaluate this choice, we compare QC-MC and QC-Titan. Note, however, that QC-MC was designed to be used after a BC. Every hour, the observations considered incorrect by the two QCs are removed. The hourly mean OmBs for the PWS temperature observations are shown in Fig. 4a. The average OmB exhibits a systematic bias following a diurnal cycle: it varies from about 0.7 °C during the day to 1.3 °C during the night. The systematic warm bias could be due to site effects as the PWSs are placed close to heating environments such as windows or balconies (Bell et al., 2015; Sgoff et al., 2022). For relative humidity (Fig. 5b), the raw PWS observations exhibit a diurnal bias composed of a drier nighttime and a wetter daytime. Its variations are close to what was reported by Sgoff et al. (2022). Regarding the observations after the two QCs, the mean OmB is reduced throughout the day when using QC-MC, while it still shows a diurnal bias with QC-Titan (Fig. 4a). QC-MC takes into account a reference, which is here the model background, allowing the QC to remove spatially incoherent observation in relation to the model background. A map of quality-controlled temperature observations near Paris at 12:00 LT is shown in Fig. 4b–e. Without QC, no regional pattern could be seen in the temperature observations from PWSs, explaining a high rms deviation from StatIC PWS observations of 1.5 °C (Fig. 4b). QC-Titan has the ability to remove observations which are too far from their neighbours. After QC-Titan, every single PWS observation is closer to its neighbours, but PWS observations have a mean deviation from StatIC PWS observations of 2.0 °C. When using QC-MC, more stations are removed, but only the observations close to the background are kept (QC-MC PWS OmB threshold is equal to ±1 °C at the time of the figure). QC-MC will therefore be chosen for the next part of the study. However, QC-MC does not ensure by itself a zero bias of PWS observations: Fig. 4a shows a slight positive bias, and a BC is necessary to remove it.

Figure 4 (a) Distributions of PWS OmB for each hour of the study period: thick lines indicate the mean, thin lines indicate the mean ± 1 standard deviation. The area between thick and thin lines is shaded. (b–e) Screen-level temperature observations around Paris on 7 September 2021 at 10:00 UTC from (coloured triangles with black contours) SWSs and (coloured stars) StatIC PWSs. Coloured circles are either (b–d) PWS screen-level temperature observations or (e) T_2 m AROME background at the PWS location. PWS observations are (b) raw or quality-controlled by (c) QC-Titan, (d) QC-MC. Boxes indicate the mean deviation (MD) and the rms deviation (RMSD) of PWS observations compared with independent StatIC PWS observations.

Figure 5 Mean diurnal cycle of OmBs during the 1-month study period of (a) temperature and (b) relative humidity for (blue) raw PWS observations; for PWS observations quality-controlled with QC-MC and bias-corrected with (brown) BC-D, (green) BC-M, (orange) BC-S, (purple, only in a) BC-R; and for (black) SWS observations. Lines and colour shades are like in Fig. 4a.

To remove the remaining bias, the four BC methods described in Sect. 2.4.1 are tested, in association with the selected QC-MC. Figure 5 shows the diurnal cycle of mean OmBs for both the PWS observations after applying the four BC methods and the SWS observations. For temperature observations (Fig. 5a), mean OmBs for all choices of BC are close to zero and are on the order of magnitude of SWS observation OmBs. By construction, BC-M mean OmBs are near zero throughout the diurnal cycle: BC-R means OmBs are very close, and the addition of hourly global solar irradiance as a predictor does not seem to make any substantial improvement. For relative humidity observations (Fig. 5b), all BCs tested exhibit a reduced diurnal cycle compared to raw PWS observations. BC-D and BC-S mean OmBs are close, as they both use 24 and 48 h rolling periods to estimate biases, respectively. Still, BC-S OmBs are slightly closer to zero during the night. Once again, the BC-M mean OmBs are near 0 % throughout the diurnal cycle, like SWS OmBs.

Figure 6 As in Fig. 4b–e, around the city of Toulouse on 7 September 2021 at 21:00 UTC (23:00 LT). Coloured circles are either (a–e) PWS screen-level temperature observations or (f) T_2 m AROME background at the PWS location. PWS observations are (a) raw and quality-controlled with QC-MC, as well as bias-corrected with (b) BC-M, (c) BC-R, (d) BC-S, and (e) BC-D. Coloured squares, stars, and triangles are Toulouse Métropole PWS, StatIC PWS, and SWS screen-level temperature observations, respectively. Boxes indicate the mean deviation (MD) and the rms deviation (RMSD) of PWS observations compared with independent Toulouse Métropole PWS observations.

Figure 6 shows how the four BC methods modify PWS observations of screen-level temperature on 7 September 2021 at 23:00 LT around the city of Toulouse (centre of the map). Raw PWS temperature observations are close to observations from SWS, Toulouse Métropole PWS, and StatIC PWS networks, even if some PWS observations are warmer north-west of the city. However, BC-M and BC-R, which use monthly rolling periods to estimate biases, reduce PWS temperatures in the suburbs of the city, which is inconsistent with other observation networks. In contrast, BC-S and BC-D, which use 24 or 48 h rolling periods to estimate biases, keep observations in agreement with other observation networks: the rms deviation (RMSD) of PWS observations compared with independent Toulouse Métropole PWS observations reaches only 0.4 and 0.2 °C, respectively. T_2 m values of the AROME background (Fig. 6f) are close to Toulouse Métropole observations (RMSD of 0.3 °C), while unfortunately none are available north-east of the city to support the colder temperatures indicated by the background.

3.2 Choice of observation error covariances and thinning

In variational DA schemes, the observation error covariance matrix R has to be specified. As in D24, the part of R allocated to PWS observations is prescribed diagonally, like for the part allocated to SWS observations, and the diagonal values (σ_o) are the same as for SWS observations: 1.4 °C for T_2 m and 10 % for RH_2 m. This is done for simplicity's sake, as modelling the non-diagonal terms is an open research question (Guillet et al., 2019), and a study of sensitivity to a change in σ_o is beyond the scope of this article. The diagonal assumption demands verifying that observation errors are not correlated, and if they are, the reduction in the spatial density of the observations may diminish the correlation between observation errors.

Figure 7 Spatial Desroziers diagnostic of PWS observations of (blue) temperature and (red) relative humidity after BC-M and QC-MC on 6 September 2021 at (solid lines) 03:00 UTC and (dashed lines) 15:00 UTC. The grey dashed line indicates the number of pairs of equidistant stations (in 1 km bins).

Figure 7 shows the spatial Desroziers diagnostic for PWS observations after BC-M and QC-MC. This diagnostic is computed using OmA and OmB of these PWS observations assimilated without thinning in a cycled 3DEnVar DA experiment launched on 6 September 2021 at 00:00 UTC. This diagnostic decreases when the distance between the observations increases. Two observations being far from each other have less probability of having correlated observation errors, reducing the influence of their local environment. The spatial Desroziers diagnostic decreases up to distances of 1 to 6 km, depending on the variable and the hour of the day. To reduce the probability of having correlated observation errors, a thinning is applied, selecting one random observation per mesh in an approximately 8 km horizontally spaced regular Gaussian grid (8 km thinning hereafter). One of the limits of the method is the fact that two observations may have correlated errors independently of the distance, for example if they are subject to similar anomalous siting conditions, and these observations are not corrected or removed by the pre-processing.

Two pre-processing methods of PWS observations are selected for the DA experiments: P-M composed of BC-M, QC-MC, and an 8 km thinning and P-S composed of BC-S, QC-MC, and an 8 km thinning. One of the points to bear in mind when using pre-processing methods based on comparisons with the background of the model into which we want to assimilate these observations is that pre-processed PWS observations can overfit the model background (see the triple collocation method shown in Appendix B).

3.3 Overview of the experiments

All experiments are described in Table 1. The benefit from assimilating PWS observations of a variable X (T for T_2 m or RH for RH_2 m) is evaluated using the observing system experiment (OSE) framework, comparing a reference experiment (3DVar, 3DEnVar) to an experiment where PWS observations are assimilated (Pourret et al., 2022; D24). To explore the role of the pre-processing (or bias correction, since this is the only part that changes), PWS observations after P-M and P-S are assimilated using the 3DEnVar DA scheme (3DEnVarX and 3DEnVarXS, respectively). To explore the role of the DA scheme, PWS observations after P-M are assimilated using the 3DVar DA scheme (3DVarX). Finally, the impact of assimilating simultaneously PWS observations after P-M at the surface using the OI DA scheme and in the atmosphere with the 3DEnVar DA scheme is tested (3DEnVarX-surf). Monitor, 3DVar, and 3DEnVar experiments are the same as in D24. Except Monitor, they are cycled, which means that the assimilation of PWS observations at time t influences the background at time t+1 h, used to make the new analysis. During 1 month, 719 analyses are made, and 119 forecasts are launched.

Table 1 Overview of the experiments. The Monitor experiment was launched on 6 August 2021 at 00:00 UTC. All cycled experiments were launched on 6 September 2021 at 00:00 UTC. All experiments use the OI surface DA scheme. X corresponds to the variable of PWS observations assimilated in the atmospheric DA scheme only unless otherwise indicated: T for T_2 m or RH for RH_2 m.

4 Results of the assimilation experiments

The experiments are evaluated using statistics of 1-month OmB and OmF (observation minus forecast), where observations are systematically SWS observations. The relative evolution of the rms OmX (X being B or F) of an experiment (XP) with respect to another (CTRL) is given by

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }\text{rms}\mathrm{OmX}={\displaystyle \frac{\mathrm{rms}\left({\mathrm{OmX}}_{\mathrm{XP}}\right)-\mathrm{rms}\left({\mathrm{OmX}}_{\mathrm{CTRL}}\right)}{\mathrm{rms}\left({\mathrm{OmX}}_{\mathrm{CTRL}}\right)}}.\end{array}$$

As in D24, Δrms OmX is considered significant if 0 is not in the 95 % confidence interval around it. This interval is computed by bootstrap with the “percentile” method (scipy.stats.bootstrap function; Virtanen et al., 2020): the OmX_XP time series is randomly sampled with a replacement, giving a new time series of the same size. This is done 1000 times, providing 1000 time series and giving 1000 rms values forming a distribution, from which the confidence interval is estimated.

4.1 Impact of the bias correction with the 3DEnVar DA scheme

When PWS T_2 m observations are assimilated in 3DEnVarT or 3DEnVarTS (Table 2), T_2 m rms OmBs decrease significantly for both experiments in comparison with 3DEnVar, which is an improvement. The improvement is slightly larger with 3DEnVarT (Δrms OmB of −0.7 %) than with 3DEnVarTS (−0.6 %) over the AROME domain. The improvement is larger (−0.9 % and −0.8 %, respectively) over France, as it is the area where PWS observations are assimilated. Surface pressure rms OmB significantly increases over France for 3DEnVarT (+1.6 %) and 3DEnVarTS (+3.7 %). Other variables show no significant evolution.

Table 2 Δrms OmB (%) of surface pressure, T_2 m, RH_2 m, and 10 m zonal and meridional wind for SWSs over the 1-month study period. Negative values of an experiment (XP) with respect to another (CTRL) indicate that the backgrounds (1 h forecasts) of XP are closer than CTRL to SWS observations; i.e. improvement and positive values indicate degradation. Significant values are in bold.

Figure 8 Δrms OmF of (a) T_2 m and (b) RH_2 m for SWSs over France every 3 h between the 0 and 24 h forecast range over the 1-month study period. Negative values indicate improvement, and positive values indicate degradation. Points (crosses) indicate improvement (degradation) of XP with respect to CTRL at the 90 % statistical significance level.

For PWS RH_2 m observation assimilation in 3DEnVarRH or 3DEnVarRHS, RH_2 m rms OmBs also decrease significantly over the AROME domain (Δrms OmB of −1.0 % and −0.8 %, respectively). The improvement is larger over France (−1.5 % and −1.2 %, respectively). However, a slight but significant degradation of T_2 m rms OmBs is found in the 3DEnVarRHS experiment over both France and the AROME domain (0.8 % and 0.5 %, respectively). Other variables show no significant evolution.

For experiments assimilating PWS observations of either T_2 m or RH_2 m, the lowest rms OmBs are obtained with the P-M pre-processing for all variables over France.

Scores for forecasts up to 24 h are shown in Fig. 8. When PWS T_2 m observations are assimilated (3DEnVarT and 3DEnVarTS in Fig. 8a), an improvement in T_2 m rms OmFs is found up to 3 h of forecasts; beyond the 3 h forecast range, the rms OmF evolution is neutral to slightly degraded. Similar results are found in 3DEnVarRH and 3DEnVarRHS when PWS RH_2 m observations are assimilated. These results are consistent with Sgoff et al. (2022) results showing improvements vanishing or not significant beyond the 5 to 6 h forecast range when assimilating bias-corrected PWS observations only.

No significant evolution of rms OmBs or rms OmFs of other observing systems is noticed (not shown), which could be explained by the small impact of T_2 m and RH_2 m above the atmospheric boundary layer (Brousseau et al., 2014). Also, for observing systems such as radiosoundings, the low number of observations over France makes it difficult to show significance (D24).

4.2 Impact of the DA scheme

When compared to 3DVar, 3DVarT significantly degrades the T_2 m rms OmBs over France by +2.4 % (Table 2). For 3DVarRH compared to 3DVar, this degradation reaches +4.5 %. The degradation is very large in the analysis and remains significant between the 6 to 12 h forecast range for T_2 m and RH_2 m, respectively (Fig. 8). As was found by D24 for PWS surface pressure observation assimilation, the 3DVar DA scheme with its operational settings is not able to take advantage of these PWS observations with the pre-processing (including the thinning) which is selected.

When PWS T_2 m observations are concomitantly assimilated by the atmospheric and surface DA systems (3DEnVarT-surf), the improvement found in 3DEnVarT turns into a significant degradation for T_2 m rms OmBs (Table 2). The degradation is significant up to 9 h forecast range (Fig. 8). When PWS RH_2 m observations are concomitantly assimilated by the atmospheric and the surface DA system (3DEnVarRH-surf), RH_2 m rms OmBs and OmFs are very similar to 3DEnVarRH, showing an improvement. These large differences between the assimilation of PWS observations of T_2 m and RH_2 m in the surface DA scheme could be explained by the fact that the RH_2 m observation operator, contrary to the T_2 m observation operator, uses only variables from the lowest atmospheric model level and not surface variables.

5 Results from case studies

In the case of a strong near-surface gradient of a meteorological variable caused by a meso-β-scale meteorological phenomenon, the added value of a dense network of weather stations is likely to be clearly visible. Two cases with substantial temperature and relative humidity gradients in September 2021 are studied: a sea breeze case during a heatwave on 7 September with around 4 °C of screen-level temperature difference between two SWS observations 12 km apart and a fog case on 23 September with 40 % screen-level relative humidity differences between two SWS observations 25 km apart.

5.1 Sea breeze front during a heatwave on 7 September 2021

On 7 September 2021, at the synoptic scale, subtropical warm air came from North Africa ahead of a low-pressure system located over the Atlantic Ocean. Near the surface, the south-eastern to eastern winds carried the warm continental air mass over western France. Maximum daily temperature reached 33.7 °C at Nantes-Bouguenais SWS at 15:07 UTC and 32.2 °C at 15:00 UTC at Vannes-Séné SWS (Fig. 9a). It is 11.3 and 10.3 °C more than the 1991–2020 normal September maximum temperature of these SWSs, respectively. The daily maximum temperature at Nantes-Bouguenais was 5 °C above or equal to the September normal from 2 to 8 September, i.e. for more than 5 consecutive days, which is the WMO criterion for heatwaves.

Figure 9 Geographical situation of (a) the heatwave case on 7 September 2021 at 15:00 UTC with raw screen-level temperature observations and (b) the fog case on 23 September 2021 at 09:00 UTC with raw screen-level relative humidity observations. SWSs are indicated by coloured triangles with black contours, PWSs by coloured circles, and StatIC PWSs by coloured stars. Human crowdsourced observations of fog between 08:00 and 09:00 UTC are indicated by three encircled horizontal bold lines (fog symbol). MSG-3 satellite HRV channel reflectance observations equal to 24 % and 43 % are shown by bold and thin black edges, respectively. Light blue indicates water surfaces (here mainly the Atlantic Ocean).

Figure 10 Analyses of the temperature of the lowest atmospheric model level (around 5 m height a.g.l.) for (a) 3DVar, (b) 3DEnVar, (c) 3DVarT, and (d) 3DEnVarT experiments on 7 September 2021 at 15:00 UTC. Locations of the stations whose observations of screen-level temperature are assimilated are indicated by triangles for SWSs and circles for PWSs. Encircled areas A and B are referred to in the text.

Figure 11 Temperature increments (a–b, d–e) at the lowest atmospheric model level (around 5 m height a.g.l.) and (c, f) at the highest surface model level (ground at 0 m height a.g.l.) on 7 September 2021 at 15:00 UTC. These experiments are designated (a) 3DVar, (b, c) 3DEnVar, (d) 3DVarT, (e) 3DEnVarT, and (f) 3DEnVarT-surf experiments. Their only difference from the experiments referred to by these four names is that they all start from the same 3DVar background at 14:00 UTC. The locations of stations whose observations of screen-level temperature are assimilated are indicated by triangles for SWSs and circles for PWSs; their colours indicate temperature OmBs.

Figure 10 shows the temperature analysis at around 5 m height a.g.l. for the different experiments and both SWS and PWS screen-level temperature observations. If the height between observations and this analysis field is different, it is shown here because it is the lowest atmospheric level to which DA increments are applied. The experiments are cycled from 6 September 2021 at 00:00 UTC onwards, which accounts for the significant discrepancy observed between the analyses. The two types of pre-processing show similar results for this case (not shown).

In the afternoon, a sea breeze rises up due to the difference in temperature between the ocean and the land, cooling the coast. This cooling is observed by three SWSs and is observed well by PWSs, particularly after pre-processing (Figs. 9a and 10). Analyses from the four DA experiments (Fig. 10) are able to reproduce the coastal cooling caused by the breeze but with differences in temperature on either side of the breeze front and differences in the location of this front.

In the 3DVar experiment (Fig. 10a) the temperature analysed is very close to SWS observations of temperature. When PWS temperature observations are assimilated in 3DVarT (Fig. 10c), the temperature analysed in area A is lower than in 3DVar, which increases the deviation from the analysis of the Vannes-Séné SWS observation. In the 3DEnVar experiment (Fig. 10b), still in area A near Vannes-Séné SWS, the breeze front is located further inland, which is not in agreement with this SWS observation or PWS observations shown in Fig 10c. Also, in area B, temperature analysed in 3DEnVar is higher than SWS temperature observations. In this case, 3DEnVarT analysis seems closer than the other analyses to both PWS and SWS observations, in particular in areas A and B (Fig. 10d).

To illustrate the shape of temperature increments given by each DA scheme depending on the observations assimilated (from SWSs only or from both SWSs and PWSs), at the lowest atmospheric level for atmospheric DA schemes and at ground level for surface schemes, experiments starting from the same 3DVar operational background at 14:00 UTC are shown in Fig. 11.

Increments at the lowest atmospheric model level have lower horizontal gradients in 3DVar (Fig. 11a) with respect to 3DEnVar (Fig. 11b) and in 3DVarT (Fig. 11d) with respect to 3DEnVarT (Fig. 11e). Because the information from observations propagates at a longer distance with the settings of the 3DVar DA scheme than the 3DEnVar DA scheme (as shown in Sect. 2.1.1), increments of the same sign as the OmBs of coastal observations propagate over longer distances at sea in 3DVar and 3DVarT than in 3DEnVar or 3DEnVarT.

The assimilation of PWS observations in 3DVarT and 3DEnVarT modifies the shape and sign of the increments when compared to 3DVar and 3DEnVar, extending positive increments southwards in particular.

In area B and generally inland, increments are much lower than the neighbouring SWS OmBs in 3DEnVar, whereas it is not the case in 3DVar. As shown in Sect. 4, the 3DEnVar DA scheme gives less relative weight to surface temperature observations, which explains why observations have little impact on the temperature analysis. Assimilating PWS observations (Fig. 11e) increases the amplitude of increments, in particular inland, in 3DVarT with respect to 3DVar or in 3DEnVarT with respect to 3DEnVar.

At ground level, the shape of the ground temperature increments with the OI in the 3DEnVar experiment (Fig. 11c) is similar to the shape of the atmospheric increments in the 3DVar experiment. When assimilating PWS temperature observations with the OI, the sign of increments at the surface is generally the same as the sign of increments added in the atmosphere in 3DVarT or 3DEnVarT, and their values are closer to 3DVarT than to 3DEnVarT.

5.2 Fog on 23 September 2021

The Saône valley, in the north-west of France (Fig. 2), is regularly affected by fog during the autumn and winter seasons. Situated between the Morvan and the Jura mountainous massif, it is crossed by a few rivers, including the Saône (Fig. 9b), providing humidity near the surface. On 23 September 2021, high-pressure conditions are associated with stable air masses over France (not shown). At 09:00 UTC (Fig. 9b), MSG-3 satellite HRV observations indicate the presence of low clouds. These clouds reach the ground, causing mist and fog: SWS and PWS relative humidity observations reach 90 % to 100 %; three human crowdsourced observations indicate fog between 08:00 and 09:00 UTC; and a forward scatter sensor at Tavaux SWS, south-west of Dole, observes a meteorological optical range of 1020 m at 09:00 UTC between mist and fog. Out of the fog, between 4.5 and 5° E, relative humidity observations range from 60 % to 90 %.

Figure 12 Relative humidity at the lowest atmospheric model level (around 5 m height a.g.l.) for (a–d) 3DEnVar and (e–h) 3DEnVarRHS experiments on 23 September 2021. (a, c, e, g) Backgrounds and (b, d, f, h) analyses are shown at two consecutive assimilation times: (a, b, e, f) 09:00 UTC and (c, d, g, h) 10:00 UTC. Assimilated observations of relative humidity from (triangles) SWSs and (circles) PWSs are indicated. MSG-3 satellite HRV channel reflectance observations equal to 24 % and 43 % are shown by bold and thin black edges, respectively.

Figure 12 shows the evolution of relative humidity in the backgrounds and analyses of two experiments at two assimilation times: 09:00 and 10:00 UTC.

At 09:00 UTC, at the lowest atmospheric level, relative humidity backgrounds of 3DEnVar (Fig. 12a) and 3DEnVarRHS (Fig. 12e) experiments are different due to the cycled assimilation. High values of relative humidity simulated by 3DEnVarRHS appear closer to the MSG-3 observation of clouds, in its east–west extension, than 3DEnVar. The analysis at 09:00 UTC (Fig. 12b, f) slightly increases the relative humidity under the observed fog but is not able to substantially modify the shape of the wet area in the background in only one assimilation time step, in particular in 3DEnVar where only SWS observations are assimilated.

At 10:00 UTC, the 1 h forecast drastically reduces the wet area for both experiments (Fig. 12c, g). Whereas in 3DEnVar the area with relative humidity above 91 % is still large in its east–west extension, in 3DEnVarRHS the wettest area is located where clouds are still observed by satellite. The assimilation at 10:00 UTC increases the relative humidity in this area, for both experiments, which exceeds 94 %, as indicated by one SWS, only in 3DEnVarRHS.

This example shows the importance of rapid update cycles in DA, as one assimilation may have a low impact, but several are gradually changing the simulated relative humidity. The model plays a preponderant role in reducing the simulated fog from 09:00 to 10:00 UTC. Here this situation is illustrated with 3DEnVarRHS because this experiment is closer to the observations than 3DEnVarRH. Regarding the results of 3DEnVarRH-surf, in this case, the assimilation of PWS observations in the OI scheme at the surface does not change the shape and the intensity of the simulated fog (not shown).

6 Discussion and conclusion

This study has explored the impact of the assimilation of temperature and relative humidity observations from PWSs in the current operational AROME system. Before the traditional screening stage of a DA system, two pre-processing methods of PWS observations were designed and selected. This was done to obtain pre-processed PWS observations with statistics of differences from the AROME background close to SWS observations currently assimilated in the AROME DA system.

• Of the four bias correction methods evaluated, two have been selected, including the method of Sgoff et al. (2022) designed to reduce the diurnal bias of PWSs. Because these methods are based on statistics of differences between AROME background and observations, there is a risk of transferring model bias to PWS observations, hence the importance of assimilating anchor near-surface observations.

• Two quality control methods have been evaluated. The Titan-QC original design with three checks has been found to be detrimental when multiple PWSs exhibit biases with respect to SWS observations; however, all checks described by Båserud et al. (2020) such as the “first guess” test have not been tested. A QC adapted from Mandement and Caumont (2020) has been used to remove PWS observations judged to be erroneous based on their departures from the AROME background.

• A thinning has been applied to reduce the probability of having correlated observation errors by taking one random observation per mesh in an approximately 8 km horizontally spaced regular Gaussian grid.

Then, the OSE framework was used: the added value of PWS observation assimilation with two selected pre-processing methods and existing DA schemes has been quantified. The assimilation of pre-processed PWS temperature observations with the 3DEnVar DA scheme reduced between 0.7 and 0.9 % the rms SWS temperature OmBs in France, depending on the pre-processing method. Departures between SWS observations and forecasts are reduced up to the 3 to 6 h range. Results are similar for PWS relative humidity observations, with a reduction of 1.2 % to 1.5 % of the rms SWS relative humidity OmBs in France. These findings are in agreement with benefits found by Sgoff et al. (2022) in experiments where PWS and SWS observations were not jointly assimilated.

Finally, the assimilation of pre-processed PWS observations with the 3DEnVar DA scheme has been subjectively found to better represent the fine-scale features of a sea breeze during a heatwave event, as well as the evolution of a fog event. In these case studies, the relative weight given to near-surface observations is found to be lower in 3DEnVar compared to 3DVar, which is compensated for by the high number of PWS observations assimilated.

Two distinct bias corrections have been evaluated. The monthly bias correction (BC-M) demonstrates superior performance in reducing the departure between SWS observations and the model background. In contrast, the 2 d temporal window bias correction (BC-S) exhibits adaptive correction capabilities that have been shown to be relevant in case studies. Further tests involving different seasons would be beneficial in order to validate the results.

This study has some limitations. The OmBs used for P-M and P-S pre-processing for all experiments are computed from a monitoring experiment with a 3DVar DA scheme (Monitor). A more rigorous calculation would require a monitoring experiment to be carried out for each assimilation scheme tested. To determine the thinning distance, the spatial Desroziers diagnostic (Fig. 7) is computed from an experiment with a 3DEnVar DA scheme. This thinning distance is applied to the experiments with a 3DVar DA scheme as well. Perhaps the scores of the 3DVar experiments are degraded by an inappropriate choice of the thinning length, which could be increased in future studies. Furthermore, no testing has been specifically conducted in complex topographical areas, such as mountainous regions, where the pre-processing may deteriorate the new observations or reject them due to temperature model bias (Gouttevin et al., 2023), or in urban areas. In DA, the observation error includes two contributions: a measurement error and a representativity error. The assumption that observational errors are uncorrelated becomes less and less valid as the average spacing between observations decreases, requiring the use of thinning techniques. Microclimatic variations, for example in urban environments, are currently considered as representativity errors if we cannot model them and are therefore also considered to be representativity errors by QC algorithms based on OmB thresholds. Advances in surface modelling, the refinement of horizontal grids, and work on modelling the correlation of observation errors should gradually reduce the proportion of microclimatic variations considered as a source of error.

The joint assimilation of temperature, relative humidity, and surface pressure observations from PWSs in AROME could now be tested, as these variables describe weather structures more coherently together than separately. Efforts should be made to evaluate such a joint assimilation over a wider variety of cases and perhaps over longer study periods. In these joint experiments, the correlation of observation errors, especially between T_2 m and RH_2 m observations, should be examined. Correlated observation error issues could be addressed by modelling a non-diagonal observation error covariance matrix or by inflating observation errors. Finally, further work remains to determine the potential for assimilating observations from PWS anemometers and rain gauges, with specific pre-processing.

For operational use, beyond the recommendations of D24 which also apply in this case, this study has not attempted to modify the settings of the assimilation systems used (3DVar, 3DEnVar, or OI), and the few new settings have been chosen identically to the operational ones. However, studies of the sensitivity of these settings would be desirable given their importance in the results obtained and because they are set to obtain the best scores with the observations currently assimilated. In the atmosphere, while the 3DEnVar system is currently only being tested, there are plans to replace it rapidly with a 4DEnVar system: it would allow more observations to be included in the analysis by making it possible to assimilate time series of observations. At the surface, the replacement of the OI by a 2DEnVar system is currently under study. This would allow more observations to be assimilated at the surface and ensure greater consistency between atmospheric and surface schemes; it is also a first step before possibly coupling surface and atmospheric DA systems.

Appendix A: Relative humidity computation

The relative humidity (RH) is calculated in AROME as

$$\begin{array}{}\text{(A1)}& \mathrm{RH}={\displaystyle \frac{\mathrm{100}Pq}{e_{\mathrm{sat}}\left(\frac{R_{\mathrm{dry}}}{R_{\mathrm{vap}}}+q\left(\mathrm{1}-\frac{R_{\mathrm{dry}}}{R_{\mathrm{vap}}}\right)\right)}}.\end{array}$$

P is the pressure (in Pa); q is the specific humidity (dimensionless); and R_dry = 287.0597 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ and R_vap= 461.5250 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ are the gas constants for dry air and water vapour, respectively.

The saturation vapour pressure e_sat (in Pa) is equal to

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{e}_{\mathrm{sat}}& =\exp ({\mathit{\alpha }}_{\mathrm{l}}+({\mathit{\alpha }}_{\mathrm{s}}-{\mathit{\alpha }}_{\mathrm{l}})\mathit{\delta }-{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{l}}+({\mathit{\beta }}_{\mathrm{s}}-{\mathit{\beta }}_{\mathrm{l}})\mathit{\delta }}{T}}\\ & -({\mathit{\gamma }}_{\mathrm{l}}+({\mathit{\gamma }}_{\mathrm{s}}-{\mathit{\gamma }}_{\mathrm{l}}\left)\mathit{\delta }\right)\ln \left(T\right)),\end{array}\end{array}$$

with δ=0 if T≥273.16 K and δ=1 if T<273.16 K and

$$\begin{array}{}\text{(A3)}& {\displaystyle }\begin{array}{rl}{\mathit{\alpha }}_{\mathrm{l}}& =\ln \left(e_{\mathrm{sat}}\right(T_t\left)\right)+{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{l}}}{T_t}}+{\mathit{\gamma }}_{\mathrm{l}}\ln \left(T_t\right),{\mathit{\beta }}_{\mathrm{l}}={\displaystyle \frac{L_{\mathrm{v}}\left(T_t\right)}{R_v}}+{\mathit{\gamma }}_{\mathrm{l}}{\mathrm{T}}_{\mathrm{t}},\\ & {\mathit{\gamma }}_{\mathrm{l}}={\displaystyle \frac{C_{\mathrm{l}}-C_{\mathrm{pv}}}{R_v}},\end{array}\text{(A4)}& {\displaystyle }\begin{array}{rl}{\mathit{\alpha }}_{\mathrm{s}}& =\ln \left(e_{\mathrm{sat}}\right(T_t\left)\right)+{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{s}}}{T_t}}+{\mathit{\gamma }}_{\mathrm{s}}\ln \left(T_t\right),{\mathit{\beta }}_{\mathrm{s}}={\displaystyle \frac{L_{\mathrm{s}}\left(T_t\right)}{R_v}}+{\mathit{\gamma }}_{\mathrm{s}}{\mathrm{T}}_{\mathrm{t}},\\ & {\mathit{\gamma }}_{\mathrm{s}}={\displaystyle \frac{C_{\mathrm{s}}-C_{\mathrm{pv}}}{R_v}}.\end{array}\end{array}$$

T is the temperature (in K), T_t=273.16 K is the triple-point temperature, $L_{\mathrm{v}}\left(T_t\right)=\mathrm{2.5008}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$ is the specific latent heat of vaporization, $L_{\mathrm{s}}\left(T_t\right)=\mathrm{2.8345}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$ is the specific latent heat of sublimation, e_sat(T_t)=611.14 Pa, $C_{\mathrm{l}}=\mathrm{4218}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ is the specific heat for water in its liquid phase, $C_{\mathrm{s}}=\mathrm{2106}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ is the specific heat for water in its solid phase, and C_pv=4R_vap is the specific heat at constant pressure and constant volume for water vapour.

Appendix B: Comparison of observations and AROME background by triple collocation

PWS and SWS screen-level observations are considered collocated when separated by less than 1 km horizontally, 50 m vertically, and less than 10 min apart in time. When it happens, the AROME background from the Monitor experiment (Table 1) corresponding to this SWS observation, given by the observation operator, provides a third value of the physical variable considered (T_2 m or RH_2 m). This triple collocation method (Stoffelen, 1998) is performed over the 1-month study period, giving a sample of 536 000 (80 000) triplets of collocated values for temperature (relative humidity).

Figure B1 Boxplots showing distributions of pairwise differences in (a) T_2 m and (b) RH_2 m in triplets of collocated values (SWS observations, PWS observations, AROME background from Monitor), depending on the pre-processing (raw, P-M without thinning, P-S without thinning) applied to PWS observations.

For each variable, three triplets with three different pre-processing methods for PWS observations are tested: (i) raw, i.e. without pre-processing; (ii) P-M without thinning; and (iii) P-S without thinning. We will refer to them simply as P-M and P-S in this paragraph only. Figure B1 shows distributions of pairwise differences between the values in each triplet.

Regarding the distributions of temperature differences between PWS and SWS collocated observations, their standard deviation is substantially reduced when pre-processing is applied, from 2.3 °C (raw) to 1.3 °C with P-M or 1.3 °C with P-S. Also, mean differences are closer to 0, and the interquartile range is reduced.

Regarding the distributions of temperature differences between the pre-processed PWS observations, with both P-M and P-S, and the AROME background, the standard deviation equals 1.0 °C, which is smaller than 1.3 °C, the standard deviation of differences between the pre-processed PWS observations and SWS observations. This could show that the selected pre-processing methods of PWS observations overfit the AROME background.

Similar conclusions can be drawn for relative humidity, although the effects of PWS pre-processing do not reduce interquartile ranges of differences from SWS observations or the AROME background by the same factor as for temperature.