Recovery of the three-dimensional wind and sonic temperature data from a physically deformed sonic anemometer

A sonic anemometer reports three-dimensional (3-D) wind and sonic temperature (T_s) by measuring the time of ultrasonic signals transmitting along each of its three sonic paths, whose geometry of lengths and angles in the anemometer coordinate system was precisely determined through production calibrations and the geometry data were embedded into the sonic anemometer operating system (OS) for internal computations. If this geometry is deformed, although correctly measuring the time, the sonic anemometer continues to use its embedded geometry data for internal computations, resulting in incorrect output of 3-D wind and T_s data. However, if the geometry is remeasured (i.e., recalibrated) and to update the OS, the sonic anemometer can resume outputting correct data. In some cases, where immediate recalibration is not possible, a deformed sonic anemometer can be used because the ultrasonic signal-transmitting time is still correctly measured and the correct time can be used to recover the data through post processing. For example, in 2015, a sonic anemometer was geometrically deformed during transportation to Antarctica. Immediate deployment was critical, so the deformed sonic anemometer was used until a replacement arrived in 2016. Equations and algorithms were developed and implemented into the post-processing software to recover wind data with and without transducer-shadow correction and T_s data with crosswind correction. Post-processing used two geometric datasets, production calibration and recalibration, to recover the wind and T_s data from May 2015 to January 2016. The recovery reduced the difference of 9.60 to 8.93 ^∘C between measured and calculated T_s to 0.81 to −0.45 ^∘C, which is within the expected range, due to normal measurement errors. The recovered data were further processed to derive fluxes. As data reacquisition is time-consuming and expensive, this data-recovery approach is a cost-effective and time-saving option for similar cases. The equation development can be a reference for related topics.

1 Introduction

The three-dimensional (3-D) sonic anemometer is commonly used for both micrometeorological research and applied meteorology (Horst et al., 2015). It directly measures boundary-layer flows at high measurement rates (10 to 50 Hz) and outputs wind speeds expressed in the 3-D right-handed orthogonal anemometer coordinate system relative to its structure frame (see Appendix A, hereafter, referred as 3-D anemometer coordinate system) and sonic temperature calculated from the speed of sound (Hanafusa et al., 1982). Its outputs are commonly used to estimate the fluxes of momentum and sonic temperature and, when combined with fast-response scalar sensors, the fluxes of CO₂, H₂O, and other atmospheric constituents.

Figure 1 Diagram of the IRGASON for the three sonic measurement paths (red dash lines) along which ultrasonic signals transmit, and the three dimensional (3-D) right-handed orthogonal anemometer coordinate system (blue lines) in which 3-D wind is expressed (i.e., u₁, u₂, and u₃ are the flow speeds along the first, second, and third sonic paths, respectively. These three flow speeds are expressed as u_x, u_y, and u_z in this 3-D anemometer coordinate system). d₃ is the third sonic path length, c₃ is the measured speed of sound along the third sonic path, and U_T is the total flow vector whose magnitude is equal to $\sqrt{u_{\mathrm{3}}^{\mathrm{2}}+u_{\perp \mathrm{3}}^{\mathrm{2}}}$ or $\sqrt{u_x^{\mathrm{2}}+u_y^{\mathrm{2}}+u_z^{\mathrm{2}}}$.

It has three pairs of sonic transducers forming three sonic paths (Fig. 1), each of which is between paired sonic transducers. The three paths are situated as optimized angles for wind measurements in the 3-D anemometer coordinate system, structuring the geometry of sonic anemometer. This geometry is quantitatively defined by the path lengths and path angles that are precisely measured during production calibration. A sonic anemometer measures the time of ultrasonic signals transmitting along each path (hereafter, referred as transmitting time). In reference to the sonic path length, the transmitting time is used to calculate the speeds of flow and sound along the path, which will be detailed in Sect. 4 as follows. According to the angles of three sonic paths, the speeds from the three paths are expressed in the 3-D anemometer coordinate system for wind and as sonic temperature for air heat property.

A sonic anemometer has geometry information embedded into its operating system (OS) for internal data processing (see Appendix A), allowing output of 3-D wind and sonic temperature. However, if it is geometrically deformed from the manufacturer's setting at millimeter scales, or even smaller, due to an unexpected physical impact in transportation, installation, or other handling, the geometry embedded in the OS is not representative of the current geometry of this sonic anemometer. As a result, the anemometer no longer outputs correct wind speeds and sonic temperatures because the deformation in geometry changes the relative spatial relationship among its six sonic transducers. If, an impact displaces a transducer relative to the others, the displacement must change at least one of the sonic path lengths and one of the sonic path angles. Fortunately, if geometrical deformation is the only problem, rather than physical damage to the transducers, the sonic anemometer can, according to its working physics (Schotland, 1955), correctly perform its transmitting-time measurements. Due to the change in a sonic path length, the speeds of air flow and sound along the path are incorrectly computed because the sonic path length embedded in the OS does not match the true length when the transmitting time was measured. As a result, the incorrect speeds along with the change in any sonic path angle might cause all 3-D wind speeds as well as sonic temperature outputs to be incorrect. These incorrect outputs are recoverable because the transmitting time was correctly measured and the deformed geometry can be remeasured (i.e., recalibrated) by the manufacturer to whom the anemometer can be shipped back with care. However, the equations and algorithms for the recovery are needed if a sonic anemometer is found to be geometrically deformed in a remote site where its use has to be continued. From such a site, it could take months, seasons, or even longer for a deformed anemometer to be transported back to the manufacturer for geometry remeasurements, recalibration, and shipped back to the site. In this case, if the measurements were not continued, a measurement season or year could be easily missed.

This study demonstrates data recovery from such a case when a sonic anemometer as a component of the IRGASON (integrated CO₂/H₂O open-path gas analyzer and 3-D sonic anemometer, Campbell Scientific Inc., 2018) was geometrically deformed during transportation to the Antarctic Zhongshan Station from China in early 2015 and had to be used until its replacement arrived at the site early the next year. If the deformed sonic anemometer was not used, one measurement-year would have been missed because the only transportation of R/V Xue Long (i.e., Snow Dragon in English) from China to the Zhongshan Station served a round-trip to the site on an annual basis. More importantly, the 2015 data were also needed by related projects for collaborations. Therefore, the geometrically deformed sonic anemometer was used to acquire the 2015 data. In early 2016, the deformed anemometer was shipped, with a pair of buffer bumpers for protection, to the manufacturer of Campbell Scientific Inc. in the US for remeasurements of its geometry to update its OS (i.e., recalibration).

Using the measurements of sonic path lengths and sonic path angles for this sonic anemometer from production calibration in April 2014 before its transportation and from recalibration in March 2016 after the field use in the Zhongshan Station, this study aims to develop and verify the equations and algorithms to recover the 2015 data measured using this geometrically deformed sonic anemometer to data as if measured with the this anemometer after recalibration although actually measured before the recalibration, providing a reference to similar cases and/or related topics.

2 Site, instrumentation, and data

The observation site was in the coastal landfast sea ice area of the Zhongshan Station (69^∘22^′ S and 76^∘22^′ E), East Antarctica (Yang et al., 2016; Yu et al., 2017; Zhao et al., 2017). In this area, as influenced by the unique solar cycles, the climate is characterized by the polar night from late March to mid-July and the polar day from mid-November to January. The polar day and the polar night are inhabitable to human life, but drive atmospheric dynamics in a way that is of interest to human beings (Valkonen et al., 2008); therefore, this region has attracted scientists to measure its surface heat balance; However, these measurements are not an easy task in terms of financial support, technical infrastructure, and administrative management. As such, only a few studies on such measurements have been conducted in this region (e.g., Vihma et al., 2009; Liu et al., 2017).

The fluxes of CO₂/H₂O, heat, radiation, momentum, and atmospheric variables were measured so that the sea ice and snow surface energy budget during both melting and frozen periods can be quantified. For these measurements, the project established two open-path eddy-covariance (OPEC) flux stations in May 2015. One station (see Fig. 2) was configured with the IRGASON (SN: 1131) for the fluxes, four-component net radiometer (model: CNR4, Kipp & Zonen, Delft, the Netherlands) for net radiation and radiation fluxes; one temperature and relative humidity probe (model: HMP155A, SN: H5140031, Vaisala, Helsinki, Finland) inside a 14-plate naturally aspirated radiation shield of model 41005 for air temperature and air relative humidity; and one infrared radiometer (model: SI-111, SN: 2962, Apogee, UT, USA) for surface temperature. In early 2016, a CSAT3B (Campbell Scientific Inc., UT, USA) was added for additional data of 3-D wind and sonic temperature. This OPEC station was also equipped with a built-in barometer (model: MPXAZ6115A, Freescale Semiconductor, TX, USA) for atmospheric pressure and a built-in 107 temperature probe (model: 100K6A1A, BetaTherm, Finland) inside a 6-plate naturally aspirated radiation shield of model 41303-5A for air temperature, the IRGASON was connected to and controlled by an EC100 electronic module (SN: 1542, OS: EC100.04.10) that, in turn, was connected to and instructed by a central CR3000 Measurement and Control Datalogger (SN: 7720, OS 25) for these sensor measurements, data processing, and data output. While receiving the data output from EC100 at 10 Hz, the CR3000 also controlled and measured slow response sensors at 0.1 Hz such as the CNR4, HMP155A, and others in support to this study. EasyFlux_CR3OP (version 1.00, Campbell Scientific Inc., 2016) was used inside CR3000. The data of 3-D wind, sonic temperature, CO₂ and H₂O amounts, atmospheric pressure, diagnosis codes for the 3-D sonic anemometer and open-path infrared gas analyzer, air temperature, and relative humidity were stored 10 records per second (i.e., 10 Hz). The data from all sensors were computed and stored by the CR3000 at every half-hour interval.

Figure 2 The eddy-covariance station located in the coastal landfast sea ice area of Antarctica Zhongshan Station (69^∘22^′ S, 76^∘22^′ E). It was configured with the IRGASON integrated CO₂/H₂O open-path gas analyzer and three-dimensional sonic anemometer, CNR4 four-component net radiometer, HMP155A air temperature and relative humidity probe, and SI-111 infrared radiometer.

3 Data check and instrument diagnosis

Immediately after the station started to run, all measured values were checked. Unfortunately, the sonic temperature from the 3-D sonic anemometer was incorrect because it was around 10 ^∘C higher than the air temperature from HMP155A or 100K6A1A. Given a H₂O density of about 1.00 g m⁻³ and air temperature about −20 ^∘C, sonic temperature should be around 0.13 ^∘C higher than air temperature (see Eq. 5 in Schotanus et al., 1983) if the sonic temperature was measured, although impossible, without an error. Further diagnosis for sonic anemometer measurements found that the sonic temperature values from the three sonic paths unexpectedly deviated around −12, 5, and −7 ^∘C, respectively, as shown by device configuration (Campbell Scientific Inc., UT, USA) connected to EC100 through a notebook computer while the station was running. Apparently, the largest absolute difference in sonic temperature among the three paths reached 17 ^∘C, although the difference from an IRGASON sonic anemometer was expected to be < 1 ^∘C. Such a large unexpected absolute difference (e.g., 17 ^∘C) among the three values from the three sonic paths might be caused by the geometrical deformation of sonic anemometer. To confirm the diagnosis, the body of the IRGASON was visually examined and painting on the knuckle of side one (i.e., first sonic path) among the top three claws was found removed as it was apparently impacted (Fig. 3). Therefore, with confidence, it was concluded that the incorrect outputs of sonic temperature were caused by the geometrical deformation of sonic anemometer while being transported to Antarctica from China. The deformation also might cause the incorrect outputs of 3-D wind. Therefore, this IRGASON should have been shipped back to the manufacturer for remeasurements of its geometry to update its OS (recalibration). However, as addressed in the introduction, the 2015 data would have been missed if it was shipped back to the manufacturer at that point. To make measurements as planned, this IRGASON continued its field duty until the next round-trip of R/V Xue Long to Antarctica from China until the end of 2015 when its replacement from the manufacturer arrived at the site.

Figure 3 Painting removed where it was apparently impacted on the knuckle of side claw (first sonic path) among the top three sonic transducer claws of the IRGASON sonic anemometer (serial no.: 1131).

In early 2016, it was replaced in the field and was shipped back to the manufacturer, where it was remeasured for sonic geometry in the recalibration process in March. The remeasurements verified our diagnosis conclusion that the IRGASON sonic anemometer was geometrically deformed (see Table A1 in Appendix A). Therefore, the 2015 data from this sonic anemometer needed to be recovered as if measured by the same anemometer after recalibration, although the data were acquired from the measurements before the recalibration.

4 Algorithm to recover the data of 3-D wind and sonic temperature

An IRGASON sonic anemometer measures wind flows along its three non-orthogonal sonic paths (i.e., the three sonic paths non-orthogonally situated in relation to each other, see Fig. 1), each of which is between a pair of sonic transducers. Sensing each other in each sonic path, the pair separately pulse two ultrasonic signals in opposite directions at the same time. The signal pulsed by the transducer facing the air flow direction along the sonic path takes less time to be sensed by its paired one than the one pulsed by the transducer against the air flow direction. In a path, the transmitting time of the ultrasonic signal upward [t_ui where subscript i can be 1, 2, or 3, denoting the sequential order of sonic path (Fig. 1). This subscript denotes the same variable throughout] and downward (t_di) are measured by the sonic anemometer (Hanafusa, 1982; Foken, 2017). In the case shown in Fig. 1 for the third sonic path, or i=3, the transmitting time of ultrasonic signal upward in the path is given by the following equation:

$$\begin{array}{}\text{(1)}& t_{u\mathrm{3}}={\displaystyle \frac{d_{\mathrm{3}}}{c_{\mathrm{3}}+u_{\mathrm{3}}}},\end{array}$$

where, along the third sonic path, d₃ is its length precisely measured during production or recalibration process using a coordinate measurement machine (CMM), c₃ is the speed of sound, and u₃ is the speed of air flow (Fig. 1); and the transmitting time of ultrasonic signal downward is given by the following equation:

$$\begin{array}{}\text{(2)}& t_{d\mathrm{3}}={\displaystyle \frac{d_{\mathrm{3}}}{c_{\mathrm{3}}-u_{\mathrm{3}}}}.\end{array}$$

4.1 Recovery of 3-D wind data

4.1.1 Algorithm of sonic anemometer to output the 3-D wind data

Equations (1) and (2) lead to

$$\begin{array}{}\text{(3)}& u_{\mathrm{3}}={\displaystyle \frac{d_{\mathrm{3}}}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{1}}{t_{u\mathrm{3}}}}-{\displaystyle \frac{\mathrm{1}}{t_{d\mathrm{3}}}}\right].\end{array}$$

Using the same procedure, u₁ and u₂ (see Fig. 1) can be derived as the same form. In reference to Eq. (3), the equation for u_i; where $i=\mathrm{1},\mathrm{2}$, or 3; can be expressed as follows:

$$\begin{array}{}\text{(4)}& u_i={\displaystyle \frac{d_i}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{u}i}}}-{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{d}i}}}\right].\end{array}$$

Similar to d₃, d₁ and d₂ are also precisely measured using CMM. The three flow speeds of u_i (i=1, 2, or 3) from the three non-orthogonal paths are expressed in the 3-D anemometer coordinate system of x, y, and z; where x and y are the horizontal coordinate axes and z is the vertical axis; and through a transform matrix A as the 3-D wind speeds (u_x, u_y, and u_z) commonly used in practical applications:

$$\begin{array}{}\text{(5)}& \left[\begin{array}{l}{u}_x\\ u_y\\ u_z\end{array}\right]=\mathbf{A}\left[\begin{array}{l}{u}_{\mathrm{1}}\\ u_{\mathrm{2}}\\ u_{\mathrm{3}}\end{array}\right],\end{array}$$

where the 3-D anemometer coordinate system (see Figs. 1 and A1) is defined by its origin at the center of sonic measurement volume, the u_x−u_y plane, parallel to the imagery plane, leveled by a built-in bubble in the anemometer structure, and the u_y−u_z plane through the first sonic path and A is a 3×3 matrix constructed using precisely measured geometry of the sonic paths in angles relative to the 3-D anemometer coordinate system (see its derivations in Appendix A). Matrix A is unique for each sonic anemometer and is embedded in its OS; therefore, the 3-D wind data outputted from the anemometer are the three components of u_x, u_y and u_z in the 3-D anemometer coordinate system.

Due to shadowing from the sonic transducer itself (transducer shadowing), the measured u_i is assumed to be lower than its true value in magnitude (Wyngaard and Zhang, 1985; Kaimal and Finnigan, 1994). As denoted by u_{Ti_n} where subscript T indicates “True” and subscript n indicates that u_{Ti_n} was estimated from n counts of iterations of transducer-shadow correction as shown in Appendix B, this true value is assumed to be approached through the transducer-shadow correction from u_i. Now, the shadow correction was implemented as an option if the OS of EC100 for the IRGASON sonic anemometer is version 5 or newer. Therefore, depending on the option, Eq. (5) alternatively can be expressed as follows:

$$\begin{array}{}\text{(6)}& \left[\begin{array}{l}{u}_x\\ u_y\\ u_z\end{array}\right]=\mathbf{A}\left[\begin{array}{l}{u}_{\mathrm{T}\mathrm{1}\mathit{\_}n}\\ u_{\mathrm{T}\mathrm{2}\mathit{\_}n}\\ u_{\mathrm{T}\mathrm{3}\mathit{\_}n}\end{array}\right].\end{array}$$

Following Host et al. (2015), based on Wyngaard and Zhang (1985), the correction equation for the sonic transducer size and sonic path geometry of the IRGASON sonic anemometer is given by

$$\begin{array}{}\text{(7)}& u_{Ti\mathit{\_}\mathrm{1}}={\displaystyle \frac{u_i}{\mathrm{0.84}+\mathrm{0.16}\sin {\mathit{\alpha }}_i}},\end{array}$$

where α_i is the angle of the total wind vector to the wind vector along sonic path i and is unknown before the two vectors are estimated, but, referencing Figs. 1 and 4, the sinα_i in Eq. (7) can be alternatively expressed as a function of flow speed values to lead Eq. (7) as follows:

$$\begin{array}{}\text{(8)}& u_{Ti}={\displaystyle \frac{u_i}{\mathrm{0.84}+\mathrm{0.16}\frac{\sqrt{U_{\mathrm{T}}^{\mathrm{2}}-u_{Ti}^{\mathrm{2}}}}{U_{\mathrm{T}}}}},\end{array}$$

where U_T is the magnitude of total true wind vector, given by

$$\begin{array}{}\text{(9)}& U_{\mathrm{T}}=\sqrt{u_x^{\mathrm{2}}+u_y^{\mathrm{2}}+u_z^{\mathrm{2}}}.\end{array}$$

In Eq. (8), all independent variables are actually related to the variables in Eq. (5). As such, using this equation, u_Ti can be computed; however, there are two inconvenient issues in this equation application to transducer-shadow corrections: (1) an analytical solution for u_Ti is not easily available because u_Ti is in a second order term under a square root in the right side of Eq. (8), although u_Ti is analytically expressed in its left side and (2) U_T is not available either because u_x, u_y, and u_z are derived from u₁, u₂, and u₃ before the transducer-shadow corrections. Fortunately, the corrections are small in magnitude, as shown in Eq. (8); therefore, u_i is closed to u_Ti. As a result, u_x, u_y, and u_z from Eq. (5) are close to those from Eq. (6). Accordingly, an iteration algorithm may be the right approach to the corrections using Eq. (8), or for the estimation of u_Ti.

Figure 4 Sonic transducer shadowing along the ith (i=1, 2, or 3) sonic path between the two sonic transducers, u_i is the measured magnitude of flow vector whose true magnitude is u_Ti; u_⊥i is the flow speed normal to the ith sonic path; u_x, u_y, and u_z are the wind speeds expressed in the three-dimensional orthogonal anemometer coordinate system; and α_i is the angle between sonic path i and the total flow vector (U_T) equal to $\sqrt{u_i^{\mathrm{2}}+u_{\perp i}^{\mathrm{2}}}$ or $\sqrt{u_x^{\mathrm{2}}+u_y^{\mathrm{2}}+u_z^{\mathrm{2}}}$. See Wyngaard and Zhang (1985) and Kaimal and Finnigan (1994) for the equation to calculate u_Ti.

For the first iteration, u_Ti in the right side of Eq. (8) could be replaced with u_i as its estimation. Given that U_T should be calculated using u_x, u_y, and u_z from Eq. (6), before the transducer-shadow corrections, U_T can be estimated using u_x, u_y, and u_z from Eq. (5); see Appendix B: Iteration algorithm for sonic transducer-shadow corrections. The iterations ensure that the difference in u_x, u_y, or u_z between the last and previous iterations are <1 mm s${}^{-\mathrm{1}}\approx \mathrm{1.96}\mathit{\sigma }$ < 1, where σ is the maximum precision (i.e., standard deviation at constant wind) among u_x, u_y, and u_z (Campbell Scientific Inc., 2018). The u_{T1_n}, u_{T2_n}, and u_{T3_n} from the last interaction are finally used for Eq. (6) to compute the 3-D wind of u_x, u_y, and u_z as sonic anemometer output.

4.1.2 Procedure to recover 3-D wind data

As addressed in Eqs. (4) to (6), a sonic anemometer measures t_ui and t_di to calculate the 3-D wind of u_x, u_y, and u_z; therefore, sonic path lengths (d_i) in Eq. (4) and transform matrix A in Eqs. (5) and (6) are embedded into the OS of sonic anemometer in the manufacture processes (see the embedded data for our study sonic anemometer in Appendix A). If the anemometer was physically deformed in transportation, installation, or other handling; the sonic path lengths and sonic path angles must be changed from what they were at the time when d_i and A were embedded into its OS; therefore, d_i in Eq. (4) and sonic path angles reflected by A in Eqs. (5) and (6) are no longer valid for this anemometer. Consequently; the output of u_x, u_y, and u_z still based on embedded d_i and A from production calibration or recalibration process are erroneous. To correct the erroneous output u_x, u_y, and u_z need to be transformed back into t_ui and t_di and be recalculated using t_ui and t_di based on the true sonic path lengths and true sonic path angles at the time when t_ui and t_di were measured in the field by the sonic anemometer physically deformed away from the manufacturer's geometrical settings before its field deployment.

For the true sonic path lengths and true sonic path angles, the IRGASON (SN: 1131) was returned to the manufacturer in the way described in Sect. 3. In the same way as in the manufacture process, the lengths and angles were remeasured using CMM. The remeasured lengths are denoted by d_Ti (i=1, 2, or 3) and the remeasured angles were used to reconstruct the transform matrix A as A_T (see Appendix A). Both d_Ti and A_T are used to update the OS of this IRGASON for future field uses and to correct u_x, u_y, u_z and T_s (sonic temperature, see Sect. 4.2) that were outputted in the field before the remeasurements. The correction procedures are different for the output of u_x, u_y, u_z with or without transducer-shadow corrections.

With transducer-shadow corrections

Transfer u_x, u_y, and u_z in the 3-D anemometer coordinate system to the flow speeds along the sonic paths after transducer-shadow corrections.

$$\begin{array}{}\text{(10)}& \left[\begin{array}{l}{u}_{\mathrm{T}\mathrm{1}\mathit{\_}n}\\ u_{\mathrm{T}\mathrm{2}\mathit{\_}n}\\ u_{\mathrm{T}\mathrm{3}\mathit{\_}n}\end{array}\right]={\mathbf{A}}^{-\mathrm{1}}\left[\begin{array}{l}{u}_x\\ u_y\\ u_z\end{array}\right]\end{array}$$

Using Eq. (B5), flow speed along the ith sonic path before transducer-shadow correction (u_i) can be expressed as follows:

$$\begin{array}{}\text{(11)}& u_i=u_{\mathrm{T}i\mathit{\_}n}\left(\mathrm{0.84}+\mathrm{0.16}{\displaystyle \frac{\sqrt{U_{\mathrm{T}}^{\mathrm{2}}-u_{Ti\mathit{\_}m}^{\mathrm{2}}}}{U_{\mathrm{T}}}}\right),\end{array}$$

where U_T can be calculated using Eq. (9) and u_{Ti_m} can be reasonably approximated using u_{Ti_n} because u_{Ti_m} and u_{Ti_n} are close enough to ensure u_x, u_y, and u_z to converge at their measurement precision (see Appendix B). Using u_i and d_i, the time term inside the square bracket in Eq. (4) can be recovered as follows:

$$\begin{array}{}\text{(12)}& \left[{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{u}i}}}-{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{d}i}}}\right]={\displaystyle \frac{\mathrm{2}{u}_i}{d_i}},\end{array}$$

Additionally, according to Eq. (4) and using d_Ti, the speed of air flow along the ith sonic path can be recalculated as u_ci:

$$\begin{array}{}\text{(13)}& u_{ci}={\displaystyle \frac{d_{Ti}}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{u}i}}}-{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{d}i}}}\right].\end{array}$$

Further replacing u_i with u_ci in the iteration algorithm for sonic transducer-shadow corrections in Appendix B, u_ci is corrected for transducer-shadowing as u_{cTi_n}. Using Eq. (6), the recovered vector of 3-D wind in the 3-D anemometer coordinate system [u_cx u_cy u_cz]^′ can be expressed as follows:

$$\begin{array}{}\text{(14)}& \left[\begin{array}{l}{u}_{cx}\\ u_{cy}\\ u_{cz}\end{array}\right]={\mathbf{A}}_{\mathrm{T}}\left[\begin{array}{l}{u}_{cT\mathrm{1}\mathit{\_}n}\\ u_{cT\mathrm{2}\mathit{\_}n}\\ u_{cT\mathrm{3}\mathit{\_}n}\end{array}\right].\end{array}$$

Without transducer-shadow corrections

Transfer u_x, u_y, and u_z in the 3-D anemometer coordinate system to the flow speeds along individual sonic paths.

$$\begin{array}{}\text{(15)}& \left[\begin{array}{l}{u}_{\mathrm{1}}\\ u_{\mathrm{2}}\\ u_{\mathrm{3}}\end{array}\right]={\mathbf{A}}^{-\mathrm{1}}\left[\begin{array}{l}{u}_x\\ u_y\\ u_z\end{array}\right]\end{array}$$

Using Eqs. (12) and (13), the speed of flow along the ith sonic path (u_ci) is recalculated (i.e., recovered). Based on Eq. (5), the recovered speeds of flow along the three sonic paths can be expressed in the 3-D anemometer coordinate system as follows:

$$\begin{array}{}\text{(16)}& \left[\begin{array}{l}{u}_{cx}\\ u_{cy}\\ u_{cz}\end{array}\right]={\mathbf{A}}_{\mathrm{T}}\left[\begin{array}{l}{u}_{c\mathrm{1}}\\ u_{c\mathrm{2}}\\ u_{c\mathrm{3}}\end{array}\right].\end{array}$$

4.2 Recover sonic temperature data

4.2.1 Algorithm of sonic anemometer to output sonic temperature

Equations (1) and (2) also lead to

$$\begin{array}{}\text{(17)}& c_{\mathrm{3}}={\displaystyle \frac{d_{\mathrm{3}}}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{1}}{t_{u\mathrm{3}}}}+{\displaystyle \frac{\mathrm{1}}{t_{d\mathrm{3}}}}\right].\end{array}$$

Using the same procedure, c₁ and c₂ (see Figs. 1 and 5) can be derived as the same form. In reference to Eq. (17), the equation for c_i, where subscript i=1, 2, or 3; can be expressed as follows:

$$\begin{array}{}\text{(18)}& c_i={\displaystyle \frac{d_i}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{u}i}}}+{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{d}i}}}\right]\end{array}$$

Here, c_i is the measured speed of sound along the sonic path i (see Fig. 5). When the crosswind (u_⊥i), or wind normal to the sonic path i, is zero; c_i is the true speed of sound (c_0i where subscript 0 indicates the speed of sound at crosswind speed equal to zero). Unfortunately, crosswind is rarely zero and c_i needs to be corrected to c_0i. According to Figs. 1 and 5, the true speed of sound is given by

$$\begin{array}{}\text{(19)}& c_{\mathrm{0}i}={\displaystyle \frac{c_i}{\cos {\mathit{\alpha }}_i}}={\displaystyle \frac{c_i}{c_i/\sqrt{c_i^{\mathrm{2}}+u_{\perp i}^{\mathrm{2}}}}}=\sqrt{c_i^{\mathrm{2}}+u_{\perp i}^{\mathrm{2}}}.\end{array}$$

Referencing the diagram for wind vectors in the left side of Fig. 5, this equation can be expressed as follows:

$$\begin{array}{}\text{(20)}& c_{\mathrm{0}i}^{\mathrm{2}}=c_i^{\mathrm{2}}+U_{\mathrm{T}}^{\mathrm{2}}-u_{Ti}^{\mathrm{2}},\end{array}$$

According to the definition of sonic temperature (Kaimal and Finnigan, 1994), the sonic temperature (K) along the ith sonic path (T_si) should be expressed as follows:

$$\begin{array}{}\text{(21)}& T_{\mathrm{s}i}={\displaystyle \frac{c_{\mathrm{0}i}^{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{d}}{R}_{\mathrm{d}}}},\end{array}$$

where γ_d (1.4003) is the ratio of dry-air-specific heat at constant pressure (1004 J K⁻¹ kg⁻¹) to dry-air-specific heat at constant volume (717 J K⁻¹ kg⁻¹) and R_d is gas constant for dry air (287.04 J K⁻¹ kg⁻¹). The sonic temperature outputted from the sonic anemometer (T_s in ^∘C) is the average from the three sonic paths (van Dijk, 2002), given by

$$\begin{array}{}\text{(22)}& T_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\sum _{i=\mathrm{1}}^{\mathrm{3}}{T}_{\mathrm{s}i}-\mathrm{273.15}={\displaystyle \frac{\mathrm{1}}{\mathrm{3}{\mathit{\gamma }}_d{R}_{\mathrm{d}}}}\sum _{i=\mathrm{1}}^{\mathrm{3}}{c}_{\mathrm{0}i}^{\mathrm{2}}-\mathrm{273.15}.\end{array}$$

Substituting c_0i with Eq. (20) and then substituting c_i with Eq. (18), T_s can be expressed as follows:

$$\begin{array}{ll}{\displaystyle }{T}_{\mathrm{s}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{3}{\mathit{\gamma }}_{\mathrm{d}}{R}_{\mathrm{d}}}}\left\{\sum _{i=\mathrm{1}}^{\mathrm{3}}\left[{\displaystyle \frac{d_i^{\mathrm{2}}}{\mathrm{4}}}{\left({\displaystyle \frac{\mathrm{1}}{t_{\mathrm{u}i}}}+{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{d}i}}}\right)}^{\mathrm{2}}-u_{Ti}^{\mathrm{2}}\right]\right.\\ \text{(23)}& {\displaystyle }& {\displaystyle }\left.+\mathrm{3}{U}_{\mathrm{T}}^{\mathrm{2}}\right\}-\mathrm{273.15}.\end{array}$$

Figure 5 Crosswind on the speed of sound. Along the ith (i=1, 2, or 3) sonic path between the two sonic transducers, u_i is the measured magnitude of flow vector whose true magnitude is u_Ti, and c_i is measured speed of sound; u_⊥i is the crosswind vector normal to sonic path i; U_T is the magnitude of total flow vector whose magnitude is equal to $\sqrt{u_i^{\mathrm{2}}+u_{\perp i}^{\mathrm{2}}}$ or $\sqrt{u_x^{\mathrm{2}}+u_y^{\mathrm{2}}+u_z^{\mathrm{2}}}$, where u_x, u_y, and u_z are the wind speeds in the three-dimensional right-handed orthogonal anemometer coordinate systems; c_0i is the speed of sound at crosswind equal to zero; and α_i is the angle between sonic path i and the total flow vector.

4.2.2 Procedure to recover sonic temperature data

Equation (23) indicates that, given d_i, a sonic anemometer estimates sonic temperature using its measured transmitting time of t_ui and t_di, the flow speeds along the sonic paths (u_i or u_Ti if corrected for transducer shadowing) that are also calculated from t_ui and t_di (see Eq. 4), and the resultant wind speed (U_T, i.e., the total wind) computed using Eq. (9), inside which the three wind components in the 3-D anemometer coordinate system are transformed from u_i using A, as explained by Eq. (5), without transducer-shadow corrections or from u_Ti also using A as explained by Eq. (6), with transducer-shadow corrections. As discussed in Sect. 4.1.2, when a sonic anemometer is geometrically deformed in an incident, the sonic path lengths and sonic path angles may be changed from what they were at the time when d_i and A were embedded into its OS; therefore, d_i in Eq. (23) and A in Eqs. (5) and (6) for u_i∕u_Ti and U_T in Eq. (23) are no longer valid for this sonic anemometer. As a result, its output of u_x, u_y, u_z, and T_s still based on embedded d_i and A must not be representative to the field wind and sonic temperature to be measured. In Sect. 4.1, the procedure to recover 3-D wind data was developed using remeasured sonic path lengths (d_Ti) and redetermined sonic path angles for A_T. The procedure to recover sonic temperature data also needs to be developed using d_Ti and recovered 3-D wind data in this section.

Based on Eq. (20), the recovered speed of sound from the sonic path i after crosswind corrections (c_c0i) can be expressed as follows:

$$\begin{array}{}\text{(24)}& c_{c\mathrm{0}i}^{\mathrm{2}}=c_{ci}^{\mathrm{2}}+U_{cT}^{\mathrm{2}}-u_{cTi}^{\mathrm{2}},\end{array}$$

where c_ci is the recovered speed of sound along sonic path i and $U_{cT}=\sqrt{u_{cx}^{\mathrm{2}}+u_{cy}^{\mathrm{2}}+u_{cz}^{\mathrm{2}}}$. After replacement of $c_{\mathrm{0}i}^{\mathrm{2}}$ with $c_{c\mathrm{0}i}^{\mathrm{2}}$ in Eq. (22), the recovered sonic temperature (T_cs in ^∘C) can be written as follows:

$$\begin{array}{}\text{(25)}& T_{cs}={\displaystyle \frac{\mathrm{1}}{\mathrm{3}{\mathit{\gamma }}_{\mathrm{d}}{R}_{\mathrm{d}}}}\sum _{i=\mathrm{1}}^{\mathrm{3}}{c}_{c\mathrm{0}i}^{\mathrm{2}}-\mathrm{273.15}.\end{array}$$

Now, the term of $c_{c\mathrm{0}i}^{\mathrm{2}}$ needs to be derived. Subtracting Eq. (20) from (24) leads to

$$\begin{array}{}\text{(26)}& c_{c\mathrm{0}i}^{\mathrm{2}}=c_{\mathrm{0}i}^{\mathrm{2}}+\left(c_{ci}^{\mathrm{2}}-c_i^{\mathrm{2}}\right)+\left(U_{cT}^{\mathrm{2}}-U_{\mathrm{T}}^{\mathrm{2}}\right)-\left(u_{cTi}^{\mathrm{2}}-u_{Ti}^{\mathrm{2}}\right).\end{array}$$

Using this equation to substitute $c_{c\mathrm{0}i}^{\mathrm{2}}$ in Eq. (25), denoting $U_{cT}^{\mathrm{2}}-U_{\mathrm{T}}^{\mathrm{2}}$ with $\mathrm{\Delta }{U}_{cT}^{\mathrm{2}}$ and denoting $u_{cTi}^{\mathrm{2}}-u_{Ti}^{\mathrm{2}}$ with $\mathrm{\Delta }{u}_{cTi}^{\mathrm{2}}$ leads to

$$\begin{array}{}\text{(27)}& T_{cs}=T_{\mathrm{s}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}{\mathit{\gamma }}_{\mathrm{d}}{R}_{\mathrm{d}}}}\sum _{i=\mathrm{1}}^{\mathrm{3}}\left[\left(c_{ci}^{\mathrm{2}}-c_i^{\mathrm{2}}\right)+\mathrm{\Delta }{U}_{cT}^{\mathrm{2}}-\mathrm{\Delta }{u}_{cTi}^{\mathrm{2}}\right].\end{array}$$

In this equation, the term of $c_{ci}^{\mathrm{2}}-c_i^{\mathrm{2}}$ is still unknown. Based on Eq. (18), $c_{ci}^{\mathrm{2}}$is given by

$$\begin{array}{}\text{(28)}& c_{ci}^{\mathrm{2}}={\displaystyle \frac{d_{Ti}^{\mathrm{2}}}{\mathrm{4}}}{\left[{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{u}i}}}+{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{d}i}}}\right]}^{\mathrm{2}}.\end{array}$$

Accordingly, the unknown term is given by

$$\begin{array}{ll}{\displaystyle }{c}_{ci}^{\mathrm{2}}-c_i^{\mathrm{2}}& {\displaystyle }={\displaystyle \frac{d_{Ti}^{\mathrm{2}}}{\mathrm{4}}}{\left[{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{u}i}}}+{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{d}i}}}\right]}^{\mathrm{2}}-{\displaystyle \frac{d_i^{\mathrm{2}}}{\mathrm{4}}}{\left[{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{u}i}}}+{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{d}i}}}\right]}^{\mathrm{2}}\\ \text{(29)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left[{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{u}i}}}+{\displaystyle \frac{\mathrm{1}}{t_{\mathrm{d}i}}}\right]}^{\mathrm{2}}\left(d_{Ti}^{\mathrm{2}}-d_i^{\mathrm{2}}\right)=c_i^{\mathrm{2}}{\displaystyle \frac{\mathrm{\Delta }{d}_{Ti}^{\mathrm{2}}}{d_i^{\mathrm{2}}}}.\end{array}$$

In this equation, the only unknown variable is $c_i^{\mathrm{2}}$. Based on Eq. (20), this equation can be expressed as follows:

$$\begin{array}{}\text{(30)}& c_{ci}^{\mathrm{2}}-c_i^{\mathrm{2}}=\left(c_{\mathrm{0}i}^{\mathrm{2}}-U_{\mathrm{T}}^{\mathrm{2}}+u_{Ti}^{\mathrm{2}}\right){\displaystyle \frac{\mathrm{\Delta }{d}_{Ti}^{\mathrm{2}}}{d_i^{\mathrm{2}}}}.\end{array}$$

In the right side of this equation, $c_{\mathrm{0}i}^{\mathrm{2}}$ is the only unknown. However, the whole term in the right side of Eq. (30) mathematically is a differential term in which $c_{\mathrm{0}i}^{\mathrm{2}}$can be reasonably approximated using its neighbor value, as close as possible to $c_{\mathrm{0}i}^{\mathrm{2}}$. The average of $c_{\mathrm{01}}^{\mathrm{2}},c_{\mathrm{02}}^{\mathrm{2}}$, and $c_{\mathrm{03}}^{\mathrm{2}}$ can be calculated from Eq. (22) because T_s is an output variable of the sonic anemometer. Without a measurement error and random error, the three c_0i should be the same, independent of flow speed, because they are the true speed of sound instead of measured speed of sound along an individual sonic path (Schotanus et al., 1983; Liu et al., 2001); Therefore, $c_{\mathrm{0}i}^{\mathrm{2}}$ can be reasonably approximated using the average of three $c_{\mathrm{0}i}^{\mathrm{2}}$ as $c_{\mathrm{0}}^{\mathrm{2}}$, given by

$$\begin{array}{}\text{(31)}& c_{ci}^{\mathrm{2}}-c_i^{\mathrm{2}}=\left(c_{\mathrm{0}}^{\mathrm{2}}-U_{\mathrm{T}}^{\mathrm{2}}+u_{Ti}^{\mathrm{2}}\right){\displaystyle \frac{\mathrm{\Delta }{d}_{Ti}^{\mathrm{2}}}{d_i^{\mathrm{2}}}},\end{array}$$

where $c_{\mathrm{0}}^{\mathrm{2}}$ can be computed from Eq. (22) as follows:

$$\begin{array}{}\text{(32)}& c_{\mathrm{0}}^{\mathrm{2}}={\mathit{\gamma }}_{\mathrm{d}}{R}_{\mathrm{d}}\left(T_{\mathrm{s}}+\mathrm{273.15}\right).\end{array}$$

Due to the replacement of $c_{\mathrm{0}i}^{\mathrm{2}}$ with $c_{\mathrm{0}}^{\mathrm{2}}$, the relative error of the whole term in the right side of Eq. (31) would be < 4 %, even if the variability in sonic temperature due to the difference among $c_{\mathrm{0}i}^{\mathrm{2}}$ values reaches 10 ^∘C at an air temperature of −30 ^∘C without wind (i.e., U_T=0 and u_Ti=0), which would be the worst case. Substituting the term of $c_{ci}^{\mathrm{2}}-c_i^{\mathrm{2}}$ in Eq. (27) with Eq. (31) leads to

$$\begin{array}{ll}{\displaystyle }{T}_{cs}& {\displaystyle }=T_{\mathrm{s}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}{\mathit{\gamma }}_{\mathrm{d}}{R}_{\mathrm{d}}}}\sum _{i=\mathrm{1}}^{\mathrm{3}}\left[\left(c_{\mathrm{0}}^{\mathrm{2}}-U_{\mathrm{T}}^{\mathrm{2}}+u_{Ti}^{\mathrm{2}}\right){\displaystyle \frac{\mathrm{\Delta }{d}_{Ti}^{\mathrm{2}}}{d_i^{\mathrm{2}}}}\right.\\ \text{(33)}& {\displaystyle }& {\displaystyle }\left.+\mathrm{\Delta }{U}_{cT}^{\mathrm{2}}-\mathrm{\Delta }{u}_{cTi}^{\mathrm{2}}\right].\end{array}$$

In the right side of this equation, the whole term after T_s is the sonic temperature recovery term.

5 Application

For our case without a transducer-shadow correction, Eqs. (15), (12), (13), and (16) were sequentially used to recover the 3-D wind data. In a case of transducer-shadow correction in option, Eqs. (10) to (16) are used. Based on the data of 3-D wind from the recovery process, Eqs. (9), (32), and (33) were used to recover the sonic temperature data. The whole recovery processes large data files (10 records per second), not only using these equations, but also operating the matrixes (A3) to (A5) (see Appendix A) for Eqs. (15) and (16) along with the data of sonic paths lengths in Table A1 for Eqs. (12) and (13). Apparently, the recovery process is a huge work load in computation. As such, these equations, matrixes, and data were implemented into a software package: “Sonic Data Recovery for IRGASON/CSAT3/A/B Used in Geometrical Deformation after Production/Calibration” (Appendix C and Fig. 6). As long as the path lengths and matrixes from production/calibration and from recalibration are input into the software as instructed by the interface (Appendix C), the software automatically recover the data in batches.

Figure 6 Dialogue interface of software: sonic data recovery for IRGASON/CSAT3/A/B used in geometrical deformation after production and calibration.

6 Verification

In our station, an additional anemometer for wind was not under deployment when this study IRGASON was used in its deformed state; therefore, no data were available to verify the recovered 3-D wind data. However, the algorithms as addressed using Eqs. (10) to (16) to recover the 3-D wind data are solid without any estimation and the recovered 3-D wind data are not necessary for verification.

Figure 7 Verification of sonic temperature (T_s) recovered against calculated (see Appendix D) from the air temperature (T), relative humidity (RH), and atmospheric pressure (P) that were measured using a HMP155A air temperature and relative humidity probe as well as the IRGASON built-in barometer. Blue curves: T_s measured by the IRGASON sonic anemometer in geometrical deformation (raw T_s); red curves: T_s recovered from raw T_s using Eq. (33); grey curves: T_s recovered also from raw T_s using Eq. (40) (i.e., adjusted Eq. 33); and green curves: T_s calculated from T, RH and P.

Fortunately, the data to verify sonic temperature are available in this station. Air temperature, relative humidity, and atmospheric pressure were measured using research grade sensors of the HMP155A and IRGASON built-in barometer and the data of these variables also stored at 10 Hz (10 records per second). These data can be used to estimate the sonic temperature (see Appendix D: Sonic temperature from air temperature, relative humidity, and atmospheric pressure). The recovered data of sonic temperature using Eq. (33) were compared to the calculated sonic temperature over the range of sonic temperature for three representative values: $-\mathrm{20.01}\pm \mathrm{0.14}$ ^∘C in Fig. 7a, $-\mathrm{9.06}\pm \mathrm{0.13}$ ^∘C in Fig. 7b, and $-\mathrm{1.90}\pm \mathrm{0.22}$ ^∘C in Fig. 7c. The difference between measured (i.e., unrecovered) and calculated sonic temperature values of 9.60±0.14 K in Fig. 7a, 9.53±0.17 K in Fig. 7b, and 8.93±0.24 K in Fig. 7c was narrowed to 0.99±0.14 K, 0.57±0.17 K, and $-\mathrm{0.25}\pm \mathrm{0.24}$ K, respectively, as the difference between recovered and calculated sonic temperature values. Given the accuracy of ±0.5 K in sonic temperature from the IRGASON sonic anemometer (Personal communication with Larry Jacobsen, the designer of the sonic anemometer, 2017) and the accuracy of ±0.2 ∼0.3 K in air temperature below 0 ^∘C and 1.2 % in relative humidity from HMP155A (Vaisala Corp., 2017), from which the calculated sonic temperature was derived (see Appendix D), recovered sonic temperature data can be reasonably judged as satisfactory if the difference in mean sonic temperature between recovered and calculated ranges within ±0.80 K or even wider, which could be considered a likelihood range of possible difference between correctly measured and calculated sonic temperature. As shown in Fig. 7, Eq. (33) apparently did an excellent job in recovering the sonic temperature data measured using sonic anemometer in its deformed state, but was less satisfactory in the case of Fig. 7a (i.e., 0.99±0.14 K, the difference in sonic temperature between recovered and calculated) although the range of 0.99±0.14 K was not significantly different from ±0.80 K. The less satisfactory recovery might be caused by the approximation of c_0i from c₀ that is fully valid if all c_0i are not measured by a sonic anemometer in its deformed state, but this is not the case in this study.

According to Eq. (22), it is impossible to have an individual c_0i from T_s, which is the sole output for sonic temperature from any sonic anemometer. Now, the average of $c_{\mathrm{01}}^{\mathrm{2}},c_{\mathrm{02}}^{\mathrm{2}}$, and $c_{\mathrm{03}}^{\mathrm{2}}$ is known and the changes in sonic path lengths are known. It is possible to estimate the difference among the three speeds of sound and to adjust their average ($c_{\mathrm{0}}^{\mathrm{2}})$ to $c_{\mathrm{01}}^{\mathrm{2}},c_{\mathrm{02}}^{\mathrm{2}}$, and $c_{\mathrm{03}}^{\mathrm{2}}$ in approximation, although the exact values are impossible to determine. The adjusted values can reflect the variability among $c_{\mathrm{0}i}^{\mathrm{2}}$ to some degree and are reasonably expected to improve the data recovery.

7 Adjustment

The measured speed of sound after crosswind correction (c_0i) is independent of wind speed (Schotanus et al., 1983; Liu et al., 2001) while depending on moist air density and atmospheric pressure (Barrett and Suomi, 1949). Without wind, c_0i is equal to the measured speed of sound (c_i) from sonic path i (see Eq. 19). In this case, again without wind, t_ui and t_di in Eq. (18) are the same and can be denoted by t_i. Accordingly, Eq. (18) in this case is equivalent to

$$\begin{array}{}\text{(34)}& c_{\mathrm{0}i}\equiv {\displaystyle \frac{d_i}{t_i}}.\end{array}$$

In Eq. (33), $c_{\mathrm{0}}^{\mathrm{2}}$ is the average of three squared c_0i (see Eqs. 22 and 32), but an individual c_0i is unknown; therefore, for recovery improvement, it has to be estimated from $c_{\mathrm{0}}^{\mathrm{2}}$ through a reasonable adjustment. The difference in magnitude between $c_{\mathrm{0}}^{\mathrm{2}}$ and $c_{\mathrm{0}i}^{\mathrm{2}}$ must be related to the $c_{\mathrm{0}i}^{\mathrm{2}}$ error due to the geometrical deformation of sonic anemometer. Squaring both sides of Eq. (34) leads to

$$\begin{array}{}\text{(35)}& c_{\mathrm{0}i}^{\mathrm{2}}={\displaystyle \frac{d_i^{\mathrm{2}}}{t_i^{\mathrm{2}}}}.\end{array}$$

The total differentiation of $c_{\mathrm{0}i}^{\mathrm{2}}$ is given by

$$\begin{array}{}\text{(36)}& \mathrm{\Delta }{c}_{\mathrm{0}i}^{\mathrm{2}}={\displaystyle \frac{\mathrm{2}{d}_i}{t_i^{\mathrm{2}}}}\mathrm{\Delta }{d}_i-{\displaystyle \frac{\mathrm{2}{d}_i^{\mathrm{2}}}{t_i^{\mathrm{3}}}}\mathrm{\Delta }{t}_i.\end{array}$$

Given the transmitting time is correctly measured by a sonic anemometer (i.e., Δt_i=0) even in its geometrical deformation, this equation becomes

$$\begin{array}{}\text{(37)}& \mathrm{\Delta }{c}_{\mathrm{0}i}^{\mathrm{2}}={\displaystyle \frac{\mathrm{2}{d}_i}{t_i^{\mathrm{2}}}}\mathrm{\Delta }{d}_i=c_{\mathrm{0}i}^{\mathrm{2}}{\displaystyle \frac{\mathrm{2}\mathrm{\Delta }{d}_i}{d_i}}=c_{\mathrm{0}i}^{\mathrm{2}}{\displaystyle \frac{\mathrm{2}\left(d_i-d_{Ti}\right)}{d_i}}.\end{array}$$

Mathematically in differentiation, $c_{\mathrm{0}i}^{\mathrm{2}}$ can be reasonably approximated by c₀, given by

$$\begin{array}{}\text{(38)}& \mathrm{\Delta }{c}_{\mathrm{0}i}^{\mathrm{2}}\approx \mathrm{2}{c}_{\mathrm{0}}^{\mathrm{2}}\left(\mathrm{1}-{\displaystyle \frac{d_{Ti}}{d_i}}\right)\end{array}$$

This is the error of $c_{\mathrm{0}i}^{\mathrm{2}}$ away from $c_{\mathrm{0}}^{\mathrm{2}}$. This error can be reasonably used to represent the deviation of $c_{\mathrm{0}i}^{\mathrm{2}}$ away from $c_{\mathrm{0}}^{\mathrm{2}}$. The deviations of three $c_{\mathrm{0}i}^{\mathrm{2}}$ values away from $c_{\mathrm{0}}^{\mathrm{2}}$ are the measures of variability among three $c_{\mathrm{0}i}^{\mathrm{2}}$ away from $c_{\mathrm{0}}^{\mathrm{2}}$.

Although an individual $c_{\mathrm{0}i}^{\mathrm{2}}$ is unknown, the average of three $c_{\mathrm{0}i}^{\mathrm{2}}$ is known as $c_{\mathrm{0}}^{\mathrm{2}}$. This average should be unchanged after adjustments because of the adjustment within the variability among $c_{\mathrm{0}i}^{\mathrm{2}}$ away from $c_{\mathrm{0}}^{\mathrm{2}}$. If the average of adjusted $c_{\mathrm{0}i}^{\mathrm{2}}$ is not equal to $c_{\mathrm{0}}^{\mathrm{2}}$, all adjusted $c_{\mathrm{0}i}^{\mathrm{2}}$ should be added or subtracted with the same constant to make the average of three adjusted $c_{\mathrm{0}i}^{\mathrm{2}}$ values as $c_{\mathrm{0}}^{\mathrm{2}}$, but the variability among $c_{\mathrm{0}i}^{\mathrm{2}}$ values is kept the same. This constant must be the mean of three $\mathrm{\Delta }{c}_{\mathrm{0}i}^{\mathrm{2}}$ values. Based on these analyses, the adjustment of $c_{\mathrm{0}}^{\mathrm{2}}$ to $c_{\mathrm{0}i}^{\mathrm{2}}$ can be constructed as follows:

$$\begin{array}{}\text{(39)}& c_{\mathrm{0}i}^{\mathrm{2}}\equiv c_{\mathrm{0}}^{\mathrm{2}}+\left(\mathrm{\Delta }{c}_{\mathrm{0}i}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\sum _{i=\mathrm{1}}^{\mathrm{3}}\mathrm{\Delta }{c}_{\mathrm{0}i}^{\mathrm{2}}\right).\end{array}$$

Using this equation to replace $c_{\mathrm{0}i}^{\mathrm{2}}$ in Eq. (30) and the resultant equation with this replacement then is used for $c_{ci}^{\mathrm{2}}-c_i^{\mathrm{2}}$ in Eq. (27) as follows:

$$\begin{array}{ll}{\displaystyle }{T}_{cs}& {\displaystyle }=T_{\mathrm{s}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}{\mathit{\gamma }}_{\mathrm{d}}{R}_{\mathrm{d}}}}\sum _{i=\mathrm{1}}^{\mathrm{3}}\left\{\left[c_{\mathrm{0}}^{\mathrm{2}}+\left(\mathrm{\Delta }{c}_{\mathrm{0}i}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\sum _{j=\mathrm{1}}^{\mathrm{3}}\mathrm{\Delta }{c}_{\mathrm{0}j}^{\mathrm{2}}\right)\right.\right.\\ \text{(40)}& {\displaystyle }& {\displaystyle }\left.\left.-U_{\mathrm{T}}^{\mathrm{2}}+u_{Ti}^{\mathrm{2}}\right]{\displaystyle \frac{\mathrm{\Delta }{d}_{Ti}^{\mathrm{2}}}{d_i^{\mathrm{2}}}}+\mathrm{\Delta }{U}_{cT}^{\mathrm{2}}-\mathrm{\Delta }{u}_{cTi}^{\mathrm{2}}\right\}.\end{array}$$

In the right side of this equation, the whole term after T_s is the adjusted sonic temperature recovery term.

The data recovered using Eq. (33) were recovered again using Eq. (40). Apparently, this equation did a better job than Eq. (33). The difference in sonic temperature between the recovered and calculated values was reduced to 0.81±0.14, 0.38±0.17, and $-\mathrm{0.45}\pm \mathrm{0.24}$ K, respectively, as shown from panels a to c in Fig. 7. These values for the difference fell into the range of ±0.80 K in a statistical sense. Eventually, Eq. (40) was used for data recovery and was incorporated into the software (Appendix C).

8 Discussion

8.1 Verification of 3-D wind recovery

Although not explicitly verified, the recovered 3-D wind data were implicitly verified through the verification of recovered sonic temperature data because (1) sonic temperature is more sensitive than wind speeds in ultrasonic sonic measurements (Thomas Foken, 2018, review comment for this publication) and (2) the recovery of sonic temperature data must rely on recovered 3-D wind data (Eqs. 33 and 40). According to Eqs. (3), (17), and (21), it is apparent that sonic temperature is sensitive to one order higher than wind speed to the errors in measurements of sonic path lengths and ultrasonic signal transmitting time values. If the recovered sonic temperature is within the accuracy limits of sensors, this should be realized for the wind data recovery as well (Thomas Foken 2018, review comment for this publication). Additionally, the cross wind correction for sonic temperature needs 3-D wind data (Liu et al., 2001). If 3-D wind had not been well recovered, sonic temperature data could not have been recovered satisfactorily. Therefore, the satisfactory recovery of sonic temperature data in this study implicitly verified the satisfactory recovery of 3-D wind data.

8.2 Comparability of recovered temperature to calculated sonic temperature

The recovered sonic temperature was sourced from the measurements of a fast response sonic anemometer, and the calculated sonic temperature was sourced from the measurements of a slow response air temperature and relative humidity probe as well as a barometer built into the IRGASON (see Appendix D). Therefore, the former reflected the fluctuations in the sonic temperature at high frequency, and the latter reflects the same fluctuations at lower frequency. As such, a pair of recovered and calculated sonic temperature values from simultaneous measurements (i.e., the same records in a time series data file) were not comparable. The difference between the pair is meaningless; therefore, the mean difference between recovered and calculated sonic temperature values over a half-hour period was used for their data comparison.

8.3 Recovered temperatures higher than calculated sonic temperature at lower temperatures

See Fig. 7. Compared to calculated sonic temperature, the recovered sonic temperature from Eq. (40) is 0.81±0.14 K higher at −20.01 ^∘C (Fig. 7a) and 0.38±0.17 K higher at −9.06 ^∘C (Fig. 7b); however, at −1.90 ^∘C, even $-\mathrm{0.45}\pm \mathrm{0.24}$ K lower (Fig. 7c). This trend of difference with temperature may be related to the performance of sonic anemometer at different temperature and the lower accuracy of temperature and humidity probe in a lower temperature range (Vaisala Corp., 2017).

The sonic path lengths and geometry of the sonic anemometer were measured in the manufacture environment of an air temperature around 20 ^∘C (i.e., manufacture temperature) and embedded into its OS for field applications. However, above or below the manufacture temperature, the sonic path lengths must become, due to thermo-expansion or -contraction of sonic anemometer structure, longer or shorter than those at the manufacture temperature while the length values of sonic paths inside the OS are unchanged. As a result, the sonic anemometer could under- or overestimate the speed of sound, thus sonic temperature. The under- or overestimation may be insignificant when temperature is not much above or below the manufacture temperature while the anemometer must work best around the manufacturer temperature. In this study, the working air temperature for the sonic anemometer was as low as −20 ^∘C, within which the sonic paths become shorter to some degree so that its measurement performance was possibly impacted. Although an assessment on the measurement performance of sonic anemometer at low or high air temperature could not be found in literature, overestimation of the speed of sound from a sonic anemometer at temperatures dozens of degrees below the manufacture temperature and thus sonic temperature is anticipated as shown in Fig. 7a to c.

However, at different air temperature the performance of the temperature and relative humidity probe and barometer built into the IRGASON, whose measurements are used to calculate the sonic temperature (see Appendix D), more stable than a sonic anemometer while their accuracies are the best at 20 ^∘C and become lower with temperature away from 20 ^∘C (Vaisala Corp., 2017). For example, HMP155A has an accuracy in air temperature to be ±0.1 ^∘C at 20 ^∘C and ±0.25 ^∘C at −20 ^∘C, as well as an accuracy in relative humidity (RH) of $\pm (\mathrm{1.0}+\mathrm{0.008}$RH) % at 20 ^∘C and to be $\pm (\mathrm{1.2}+\mathrm{0.012}$RH) % at −20 ^∘C. The greater disagreement between recovered and calculated sonic temperature values at lower temperature in Fig. 7a may also be due to the fact that the lower the air temperature, the lower the accuracies of HMP155A and the barometer.

8.4 Radiation on calculated sonic temperature

Compared to the recovered sonic temperature using Eq. (40), the calculated sonic temperature was 0.45±0.24 ^∘C higher over a whole period of 12:00 to 12:30 and even 0.65±0.19 ^∘C higher over a partial period of 12:15 to 12:27, which may be contributed to in part by higher incoming solar radiation of 750 W m⁻² in short-wave on the radiation shield of HMP155A (Fig. 7c). As addressed in Appendix D, the calculated sonic temperature was sourced from the measurements of air temperature and relative humidity from HMP155A, as well as atmospheric pressure from the barometer built into the IRGASON. The HMP155A housed inside a radiation shield (Fig. 2) was subject to contamination from solar radiation. A radiation shield was used to shade HMP155A from sunlight, when such a shield was used, any heat generated from the shield under sunlight and the sensor under electronic power was dissipated inefficiently (Lin et al., 2001). As a result, the air and HMP155A sensing elements inside the shield were warmer than ambient air of interest. How warm the air is inside the radiation shield depended on shield structure, ambient wind speed, and other environmental conditions (Blonquist et al., 2009). In the case of Fig. 7c at 750 W m⁻² of incoming short-wave radiation, air being a degree warmer inside the radiation shield was not unusual (Lin et al., 2001). In our study, this higher air temperature could directly cause the overestimation of calculated sonic temperature (Eq. D1 in Appendix D).

8.5 Possibility and necessity of recovering the data from a geometrically deformed sonic anemometer for fluxes

A geometrically deformed sonic anemometer outputs erroneous data. These data may be recoverable or unrecoverable, depending on the degree of deformation. If the degree is too large, the sonic anemometer cannot perform its normal measurements for the transmitting time. In this case, a Campbell sonic anemometer sets high for one to six of its first six measurement warning flags (low amplitude, high amplitude, poor signal lock, large sonic temperature difference, ultrasonic signal loss, and calibration signature error; see Table 10-2 in Campbell Scientific Inc., 2018). The geometrical deformation in sonic paths could trigger one or two flags high that indicate poor signal lock and/or ultrasonic signal loss. Regardless, in the case that any of the six warning flags from a deformed sonic anemometer were frequently, regularly, or continuously high, the erroneous data must not be recoverable (i.e., data recovery is not possible). While all six warning flags are low under normal measurement conditions, the transmitting time of ultrasonic signals along each sonic path is correctly measured and the data should be recoverable. The 3-D wind data can be recovered without uncertainty although there is little uncertainty in sonic temperature (see Eqs. 33 and 40). The subsequent question is the necessity to recover the recoverable data.

A sonic anemometer is used primarily for the fluxes of momentum and heat from the fluctuations in 3-D wind speeds and sonic temperature. If the fluctuations are not significantly influenced by the geometric deformation of the sonic anemometer, the data from this anemometer may not need recovering although the data are recoverable. The fluctuations in a wind speed component or sonic temperature are measured by variance. Therefore, this influence of sonic anemometer deformation on fluctuations in wind speed and sonic temperature can be tested through analyzing the homogeneity in variance of each wind component and sonic temperature between unrecovered and recovered data.

For this study case, the 2-day data, without a missing record or any high warning flag from 10 and 11 May 2015, were used for the analyses. After data recovery processing (Fig. 6), two datasets, unrecovered and recovered, were acquired. In the unrecovered dataset, for each wind speed component or sonic temperature, the data of 18 000 values from each half-hour were used to compute its variance ($s_k^{\mathrm{2}})$, given by

$$\begin{array}{}\text{(41)}& s_k^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{18}\mathrm{000}}}\sum _{j=\mathrm{1}}^{\mathrm{18}\mathrm{000}}{\left(x_{kj}-{\overline{x}}_k\right)}^{\mathrm{2}},\end{array}$$

where x represents u_x, u_y, u_z, or T_s; subscript j denotes the jth values in kth half-hour interval, and the upper bar indicates the average over the interval. In the recovered dataset, this variance was similarly computed and denoted by $s_{Rk}^{\mathrm{2}}$, where subscript R indicates that this variance was computed from recovered dataset. For each wind component or sonic temperature, 96 variance values were available in each datasets and 192 variance values were available in both dataset. The 192 variance values for each wind component or sonic temperature can be used to construct an F-statistic (Snedecor and Cochran, 1989) to analyze the homogeneity in variance of each wind component or sonic temperature between unrecovered and recovered data, given by

$$\begin{array}{}\text{(42)}& \sum _{k=\mathrm{1}}^{\mathrm{96}}{s}_k^{\mathrm{2}}/\sum _{k=\mathrm{1}}^{\mathrm{96}}{s}_{Rk}^{\mathrm{2}}\sim F(\mathrm{1}\mathrm{727}\mathrm{904},\mathrm{1}\mathrm{727}\mathrm{904}).\end{array}$$

From this statistic, four F values were acquired for three wind components and sonic temperature. The four F values were either > 1.00 or < 1.00, showing the inhomogeneity in variance between unrecovered and recovered data (P < 0.001), which indicates that the geometrical deformation of the sonic anemometer did significantly influence the fluctuations in each of its measured variables.

Further, using EddyPro (LI-COR Biosciences, 2016), the same datasets were used to compute two sets of sensible heat flux, latent heat flux, and CO₂ flux for each half-hour interval. One set was computed using unrecovered data and the other set from recovered data. The two sets of flux data were shown in Fig. 8. Compared to the flux from unrecovered data, the flux from recovered data was 1.5 W m⁻² lower for sensible heat (P=0.031), 0.14 W m⁻² higher for latent heat (P=0.001), and 0.08 µmol m⁻² s⁻¹ higher for CO₂ (P=0.000). These values were small in magnitude, but significant in comparison to these flux values over the ice surface in Antarctica.

Analyses of the F tests and Fig. 8 show that the data measured from a geometrically deformed sonic anemometer need to be recovered; otherwise, there were significant uncertainties in the wind speed and sonic temperature fluctuations for flux estimations.

Figure 8 Comparison of sensible heat flux, latent heat flux, and CO₂ flux from recovered data (red curves) to those from unrecovered data (blue curves). The mean difference (the green bar represents the red curve minus blue curve value) is −1.5 W m⁻² < 0 (P=0.031) for sensible heat flux, 0.14 W m⁻² > 0 (P=0.001) for latent heat flux, and 0.08 µmol m⁻² s⁻¹ > 0 (P=0.000) for CO₂ flux.

8.6 Applicability of equations and algorithms in this study

Any sonic anemometer is slender (e.g., < 1.00 cm in each diameter of six claws to hold individual sonic transducers) and as light as possible to minimize its aerodynamic resistance to air flows and to maximize its stability on supporting infrastructure (e.g., tripod) to wind momentum load, which sacrifices its durability in keeping its geometrical shape. Therefore, a sonic anemometer is easily deformed if not well cared for during transportation (e.g., the case in this study), installation, or other handling. As shown in this study, a slight geometrical deformation of sonic path length as small as millimeters or less (see Table A1 in Appendix A) could cause significant errors in 3-D wind and especially in sonic temperature. According to our recalibration experience with 3-D sonic anemometers at Campbell Scientific Inc., these cases as addressed in this study have been not unusual, but the equations and algorithms to recover the data measured by a deformed 3-D sonic anemometer were not available. As requisitions of these datasets are expensive, their recovery would be a cost-effective and time-saving option.

The equations and algorithms in this study were developed based on the measurement working physics and sonic path geometry of the IRGASON sonic anemometer. The physics is the same as those for other models of Campbell Scientific 3-D sonic anemometers in use, such as CSAT3, CSAT3A, and CSAT3B (Campbell Scientific Inc., UT, USA; Horst et al., 2015). However, the sonic path geometry of the IRGASON sonic anemometer is different from other models in the assigned azimuth angle of the first sonic path in the 3-D anemometer coordinate system. This angle was assigned as 90^∘ in the IRGASON sonic anemometer, but as 0^∘ in other models (e.g., CSAT3, CSAT3A, and CSAT3B). Even so, given the sonic path lengths and transfer matrixes of sonic anemometer that were measured and determined in the manufacture or calibration process (d_i in Eq. 12 and A in Eq. 15) and in the recalibration process after use in the geometrical deformation state (d_Ti in Eqs. 13, 33, and 40 and A_T in Eqs. 14 and 16), the equations and algorithms from this study are applicable to all models of Campbell Scientific 3-D sonic anemometers (Fig. 6) except for CSAT3 if its bugged OS version 4 is used (Burns et al., 2012). The derivation procedures and even equations based on the measurement working physics are applicable as a reference to the development of the equations and algorithms to recover the data measured using other brands of 3-D sonic anemometers that incurred deformations or to studies on similar topics.

9 Conclusion remarks

An IRGASON 3-D sonic anemometer (SN: 1131) was geometrically deformed by an impact during transportation to Antarctica from China in early 2015. To fulfill the field measurement plans for the year, it had to be deployed there in the Zhongshan Station until early 2016 when it was replaced in the field with another IRGASON provided by the manufacturer and was returned to the manufacturer, Campbell Scientific Inc., for recalibration through the remeasurement of its sonic path geometry (lengths and angles), redetermination of its transfer matrix, and an update of its operating system (OS). To recover the 3-D wind and sonic temperature data measured by this sonic anemometer in its deformed state before the recalibration, equations and algorithms were developed and implemented into a software package: “Sonic Data Recovery for IRGASON/CSAT3/A/B Used in Geometrical Deformation after Production/Calibration” (Fig. 6 and Appendix C). Given two sets of sonic path lengths and two transfer matrixes of sonic anemometer that were measured and determined in the manufacture and calibration process and also in recalibration process after the use in its deformed state, the data measured by the IRGASON 3-D sonic anemometer, even in its deformed state, were recovered as if measured by the same anemometer recalibrated immediately after its deformation.

Inside a Campbell Scientific sonic anemometer, the transducer-shadow correction for 3-D wind (Wyngaard and Zhang, 1985) is an available, programmable option for a user. However, the crosswind correction for sonic temperature (Liu et al., 2001) is internally applied as default by its OS. In a case of transducer-shadow correction in option, the 3-D wind data are recovered using Eqs. (10) to (16). If not, Eqs. (15), (12), (13), and (16) are sequentially used. Based on the data from the recovery process of 3-D wind, the sonic temperature data are recovered using Eqs. (9), (32), (38), and (40); therefore, the satisfactory recovery for both 3-D wind data and sonic temperature can be eventually reflected by the satisfactory of sonic temperature data recovery.

The software based on the equations and algorithms from this study can recover the 3-D wind data with or without the transducer-shadow correction inside the sonic anemometer and sonic temperature data with crosswind correction also inside the sonic anemometer. It was verified by comparing the recovered to calculated sonic temperature data (Appendix D). As shown in Fig. 7, the recovered data of sonic temperature using Eqs. (33) and (40) were compared to the calculated sonic temperature of three representative values over the range of measured sonic temperature from −20.01 to −1.90 ^∘C. The difference of 9.60 to 8.93 K between unrecovered and calculated sonic temperature (i.e., unrecovered minus calculated) was narrowed by Eq. (40) to 0.81 to −0.45 K (i.e., recovered minus calculated), which was satisfactory for measurements of sonic anemometer below 0 to −20 ^∘C. After verification, the software was used to recover the data measured by the IRGSON (SN: 1131) 3-D sonic anemometer in its deformed state from May 2015 to January 2016. The 8-month data were recovered using 3 days of one engineer's time. Further, using EddyPro 6.2.0 (LI-COR Inc., 2016), the recovered data were processed for the fluxes of CO₂/H₂O, sensible heat, and momentum. The data quality (Foken et al., 2012) mostly ranged in 1 to 3 and the energy closure without considering surface heat flux into ice were > 83 % when friction velocity was > 0.2 m s⁻¹. Although energy balance closure is not a good indicator for data quality (Foken et al., 2012), this closure rate is fair.

The use of a deformed 3-D sonic anemometer is a practical case. The analyses of our study case indicated that the measured fluctuations in wind speeds and sonic temperature as well as fluxes were significantly influenced by the deformation. If the data from such a use cannot be recovered, the requisition of these data is expensive and their recovery would be a cost-effective and time-saving option. The equations, algorithms, and software are applicable to all models of Campbell Scientific 3-D sonic anemometers such as CSAT3, CSAT3A, and CSAT3B that are used around the world. The derivation procedures and even equations based on the measurement working physics of sonic anemometers are applicable as a reference to the development of the equations and algorithms to manage the data measured using other brands of 3-D sonic anemometers or recover the data measured by an anemometer in its deformed state.

Appendix A: Transform matrixes

In micrometeorological applications, the wind speeds are expressed in a three-dimensional (3-D) orthogonal coordinate system of anemometer or natural wind, but a sonic anemometer measures flow velocities along its three non-orthogonal sonic paths (i.e., situated non-orthogonally from each other, see Figs. 1 and A1); therefore, for applications, the flow velocities along the three sonic paths need to be transformed into a 3-D right-handed orthogonal coordinate system in reference to the geometry of sonic anemometer, as shown in Fig. A1 (i.e., the 3-D orthogonal anemometer coordinate system). Given u_x and u_y are two horizontal velocities in the x and y direction, respectively, and u_z is vertical velocity in the z direction (Fig. A1); x, y, and z are the three coordinate axes in the 3-D orthogonal anemometer coordinate system. This system is defined with the x–y plane, parallel to the anemometer bubble-leveled plane, with the first sonic path on the y–z plane, and with origin in the center of measurement volume. A flow speed along the ith ($i=\mathrm{1},\mathrm{2}$, or 3) sonic path is a combination of component velocities of u_x, u_y, and u_z; given by

$$\begin{array}{}\text{(A1)}& {\displaystyle }{u}_i=\left(u_x\cos {\mathit{\varphi }}_i+u_y\sin {\mathit{\varphi }}_i\right)\sin {\mathit{\theta }}_i+u_z\cos {\mathit{\theta }}_i,\end{array}$$

where θ_i and ϕ_i are the zenith and azimuth angles of the ith sonic path in the 3-D orthogonal anemometer coordinate system. In this system (see Fig. A1), given the first sonic path has an azimuth angle of ϕ₁ equal to 90^∘ as fixed on the x−y plane, Eq. (A1) can be expressed in a matrix form of

$$\begin{array}{l}{\displaystyle }\left[\begin{array}{c}{u}_{\mathrm{1}}\\ u_{\mathrm{2}}\\ u_{\mathrm{3}}\end{array}\right]=\left[\begin{array}{ccc}\mathrm{0}& \sin {\mathit{\theta }}_{\mathrm{1}}& \cos {\mathit{\theta }}_{\mathrm{1}}\\ \sin {\mathit{\theta }}_{\mathrm{2}}\cos {\mathit{\varphi }}_{\mathrm{2}}& \sin {\mathit{\theta }}_{\mathrm{2}}\sin {\mathit{\varphi }}_{\mathrm{2}}& \cos {\mathit{\theta }}_{\mathrm{2}}\\ \sin {\mathit{\theta }}_{\mathrm{3}}\cos {\mathit{\varphi }}_{\mathrm{3}}& \sin {\mathit{\theta }}_{\mathrm{3}}\sin {\mathit{\varphi }}_{\mathrm{3}}& \cos {\mathit{\theta }}_{\mathrm{3}}\end{array}\right]\\ \text{(A2)}& {\displaystyle }\left[\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right]={\mathrm{A}}^{-\mathrm{1}}\left[\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right],\end{array}$$

where A is a matrix expressing the flow speeds along the three non-orthogonal sonic paths in the 3-D orthogonal anemometer coordinate system. Nominally for the sonic paths of the IRGASON, θ₁, θ₂, and θ₃ are all 30^∘ and ϕ₂ and ϕ₃ are 330 and 210^∘, respectively (see Fig. A1). Given ϕ₁=90^∘, these angles are calculated using measured data from a coordinate measurement machine and, along with the sonic path lengths, are listed in Table A1 for the IRGASON serial no. of 1131 before and after its geometrical deformation.

Figure A1 IRGASON sonic path angle geometry in the three-dimensional right-handed anemometer coordinate system of x, y, and z. Blue arrows are coordinates; a red arrow between a pair of sonic transducers is the sonic path vector whose direction is defined for air flow direction, the red arrow below the IRGASON is the projection of the corresponding sonic path vector on the x–y plane, i.e., anemometer (instrument) bubble-leveled plane. As indicated by their subscript of 1, 2, or 3 for the first, second, or third sonic path, θ₁, θ₂, and θ₃ are their zenith angles and φ₁, φ₂, and φ₃ are their azimuth angles.

Table A1 The lengths, zenith angles, and azimuth angles of sonic paths in the IRGASON (serial no.: 1131) anemometer coordinate system before and after its geometrical deformation (measured using a coordinate measurement machine on 9 September 2014 before the deformation and on 6 March 2016 after deformation).

Using the data in this table, matrix A in Eq. (A2) and its inversion A⁻¹ for this IRGASON before its geometric deformation (i.e., as used in the IRGASON OS but not valid in the field after deformation) are given as follows:

$$\begin{array}{}\text{(A3)}& {\displaystyle }\mathbf{A}=\left[\begin{array}{ccc}\mathrm{0.034785}& \mathrm{1.142665}& -\mathrm{1.183914}\\ \mathrm{1.365505}& -\mathrm{0.696580}& -\mathrm{0.660515}\\ \mathrm{0.367627}& \mathrm{0.401124}& \mathrm{0.380356}\end{array}\right],\end{array}$$

and

$$\begin{array}{}\text{(A4)}& {\displaystyle }{\mathbf{A}}^{-\mathrm{1}}=\left[\begin{array}{ccc}\mathrm{0.00000}& \mathrm{0.499023}& \mathrm{0.866589}\\ \mathrm{0.418196}& -\mathrm{0.246062}& \mathrm{0.874394}\\ -\mathrm{0.441030}& -\mathrm{0.222826}& \mathrm{0.869391}\end{array}\right].\end{array}$$

After the IRGASON geometrical deformation, matrix A became

$$\begin{array}{}\text{(A5)}& {\displaystyle }{\mathbf{A}}_{\mathrm{T}}=\left[\begin{array}{ccc}\mathrm{0.006035}& \mathrm{1.276412}& -\mathrm{1.323287}\\ \mathrm{1.363991}& -\mathrm{0.724862}& -\mathrm{0.600545}\\ \mathrm{0.368690}& \mathrm{0.417250}& \mathrm{0.345690}\end{array}\right],\end{array}$$

where subscript T indicates “True” because, after the IRGASON deformation, it should be used in the field although it was not used. The inversion of this matrix is given as follows:

$$\begin{array}{}\text{(A6)}& {\displaystyle }{\mathbf{A}}_{\mathrm{T}}^{-\mathrm{1}}=\left[\begin{array}{ccc}\mathrm{0.000000}& \mathrm{0.498879}& \mathrm{0.866672}\\ \mathrm{0.347992}& -\mathrm{0.246063}& \mathrm{0.904629}\\ -\mathrm{0.420029}& -\mathrm{0.235072}& \mathrm{0.876537}\end{array}\right].\end{array}$$

Matrixes A⁻¹, A_T, and ${\mathbf{A}}_{\mathrm{T}}^{-\mathrm{1}}$ were used for our data recovery and A was also used in the sonic anemometer OS.

Appendix B: Iteration algorithm for sonic transducer-shadow corrections

Given transform matrix A, using Eq. (5), the measured wind vector [u₁ u₂ u₃]^′ along the sonic paths is transformed to the wind vector in the 3-dimensional orthogonal anemometer coordinate system [u_x u_y u_z]^′. Subsequently, U_T is calculated using Eq. (9). Replace u_Ti with u_i under the square root in the right side of Eq. (8), an approximate equation for the first iteration is given as follows:

$$\begin{array}{}\text{(B1)}& {\displaystyle }{u}_{Ti\mathit{\_}\mathrm{1}}\approx {\displaystyle \frac{u_i}{\mathrm{0.84}+\mathrm{0.16}\frac{\sqrt{U_{\mathrm{T}}^{\mathrm{2}}-u_i^{\mathrm{2}}}}{U_{\mathrm{T}}}}},\end{array}$$

where i is 1, 2, or 3 and subscript 1 of u_Ti indicates that it is calculated from the first iteration.

First iteration

Equation (B1) is used for sonic transducer-shadow corrections in the first iteration.

Second iteration

$$\begin{array}{}\text{(B2)}& {\displaystyle }\left[\begin{array}{l}{u}_x\\ u_y\\ u_z\end{array}\right]=\mathbf{A}\left[\begin{array}{l}{u}_{T\mathrm{1}\mathit{\_}\mathrm{1}}\\ u_{T\mathrm{2}\mathit{\_}\mathrm{1}}\\ u_{T\mathrm{3}\mathit{\_}\mathrm{1}}\end{array}\right]\end{array}$$

Using Eq. (9), U_T is recalculated. Replace u_i with u_{Ti_1} under the square root in the right side of Eq. (B1), an approximate equation for the second iteration is given as follows:

$$\begin{array}{}\text{(B3)}& {\displaystyle }{u}_{Ti\mathit{\_}\mathrm{2}}={\displaystyle \frac{u_i}{\mathrm{0.84}+\mathrm{0.16}\frac{\sqrt{U_{\mathrm{T}}^{\mathrm{2}}-u_{Ti\mathit{\_}\mathrm{1}}^{\mathrm{2}}}}{U_{\mathrm{T}}}}}\end{array}$$

Third iteration

…

nth iteration

$$\begin{array}{}\text{(B4)}& {\displaystyle }\left[\begin{array}{l}{u}_x\\ u_y\\ u_z\end{array}\right]=\mathbf{A}\left[\begin{array}{l}{u}_{T\mathrm{1}\mathit{\_}m}\\ u_{T\mathrm{2}\mathit{\_}m}\\ u_{T\mathrm{3}\mathit{\_}m}\end{array}\right]\end{array}$$

where subscript $m=n-\mathrm{1}$. Using Eq. (9), U_T is also recalculated. Similar to the calculation for u_{Ti_2}, u_{Ti_n} is calculated using the following equation:

$$\begin{array}{}\text{(B5)}& {\displaystyle }{u}_{\mathrm{T}i\mathit{\_}n}={\displaystyle \frac{u_i}{\mathrm{0.84}+\mathrm{0.16}\frac{\sqrt{U_{\mathrm{T}}^{\mathrm{2}}-u_{Ti\mathit{\_}m}^{\mathrm{2}}}}{U_{\mathrm{T}}}}},\end{array}$$

to ensure that the difference in u_x, u_y, or u_z between the last and previous iterations is <1 mm s${}^{-\mathrm{1}}\approx \mathrm{1.96}\mathit{\sigma }$, where σ is the maximum precision (i.e., standard deviation at constant wind) among u_x, u_y, and u_z (Campbell Scientific Inc., 2018). Our numerical tests within the measurement ranges in u_x, u_y, and u_z concluded that the iterations mostly converged at n=2 and entirely at n≤3.

Appendix C: MATLAB code

Sonic data recovery for the IRGASON/CSAT3/A/B used in geometrical deformation after production/calibration (Code lines were formatted for readability and the electronic version of this code is available from the corresponding authors).

Note: This code can be compiled in MATLAB as an executable file: Data_recovery.exe.

% sonicdatarecovery Sonic Data Recovery for IRGASON/CSAT3/A/B Used in Geometrical Deformation after Production/Calibration

%Syntax:

function [Ux,Uy,Uz,Ts,Ts1,Ts2,Raw] = 
sonicdatarecovery(RAW)

% Inputs:

% um Measured 3-D wind speeds in the orthogonal anemometer coordinate system (OCS)

% T_s Measured sonic temperature

% A Matrix of sonic to OCS before geometrical deformation

% A_T Matrix of sonic to OCS after geometrical deformation

% di Sonic path length before geometrical deformation (i=1, 2, or 3)

% dTi Sonic path length after geometrical deformation (i=1, 2, or 3)

% Constants

shadow_correction_flag =1;  % %Shadow correction has been done (=1) or not (=0) inside OS

gama_d=1.4003; %% the ratio of dry air specific heat at constant pressure to that at constant volume

Rd =287.04; %% gas constant for dry air

RV $=\mathrm{4.61495}e-\mathrm{4}$; %% gas constant for water vapor

Av =60.064621; Bv =60.973392; Cv =60.387959;
Ah =0.000000; Bh =59.527953; Ch =63.195226;

Avt =60.074122; Bvt =64.773415; Cvt =61.227399;
Aht =0.000000; Bht =54.736084; Cht =60.766176;

% Browse to the raw data file directory to load files in a batch

hwait = waitbar(0,“Please select the file to be processed”);

pause(0.5)

[name,path] = uigetfile(“*.*”,“stabilitylect a folder”);

fname = [path name];

close(hwait);

RAW = dlmread(fname,',', 4, 1);

% Extract sonic anemometer and other meteorological data

UX = RAW(:,2); UY = RAW(:,3); UZ = RAW(:,4);

TRAW = RAW(:,5); H2O = RAW(:,8);
Temp = RAW(:,10); P = RAW(:,11);

amb_e = RV.*H2O.*(Temp +273.15);
TS_emp = (Temp +273.15).*(1+0.32*amb_e./P)−273.15;

% Load transform matrix of Eq. (A2) and data of Table A1 before geometrical deformation

The1 = ((90-Av)/180)*pi; The2 = ((90-Bv)/180)*pi;
The3 = ((90-Cv)/180)*pi;

Phi1 = ((90-Ah)/180)*pi; Phi2 = ((270+Bh)/180)*pi;
Phi3 = ((270-Ch)/180)*pi;

A_inversion = [0 sin(The1) cos(The1);
sin(The2)*cos(Phi2) sin(The2)*sin(Phi2) cos(The2);

sin(The3)*cos(Phi3) sin(The3)*sin(Phi3) cos(The3)];

A=A_inversion(−1); $d=[\mathrm{11.6486};\mathrm{11.5240};\mathrm{11.4968}]$;

% Load transform matrix of Eq. (A5) and data of Table A1 after geometrical deformation

The1 = ((90-Avt)/180)*pi; The2 = ((90-Bvt)/180)*pi;
The3 = ((90-Cvt)/180)*pi;

Phi1 = ((90-Aht)/180)*pi; Phi2 = ((270+Bht)/180)*pi;
Phi3 = ((270-Cht)/180)*pi;

AT_inversion = [0 sin(The1) cos(The1);
sin(The2)*cos(Phi2) sin(The2)*sin(Phi2) cos(The2);
sin(The3)*cos(Phi3) sin(The3)*sin(Phi3) cos(The3)];

AT = AT_inversion^∧(-1);
dT =[11.6159;11.1245;11.3548];

% Prompt data processing is in progress

hwait = waitbar(0,“Processing> > >  > > > ”)

%Recover 3-D wind data

%Get measured flow speeds along each of 3 sonic paths

[mRaw,nRaw] = size(RAW);

for i=1:mRaw;

um = [UX(i);UY(i);UZ(i)];

%With transducer-shadow corrections (TSC):

UT = (um(1)²+um(2)^∧2+um(3)^∧2)${}^{\wedge }(\mathrm{1}/\mathrm{2})$; %% Calculate the total wind magnitude

if isequal(shadow_correction_flag, 1) %% TSC has been done (=1) inside firmware

u= A_inversion*um; %% Calculate the vector of the three flow speeds using Eq. (10)

ut1(1) =u (1)/(0.84+0.16.*((UT^∧2-u
(1)^∧2)${}^{\wedge }(\mathrm{1}/\mathrm{2})$)./UT);
%% Eq. (11), recover flow speed along sonic path 1 before TSC

ut2(1) =u (2)/(0.84+0.16.*((UT^∧2-u
(2)^∧2)${}^{\wedge }(\mathrm{1}/\mathrm{2})$)./UT);
%% Eq. (11), recover flow speed along sonic path 2 before TSC

ut3(1) =u (3)/(0.84+0.16.*((UT^∧2-u
(3)^∧2)${}^{\wedge }(\mathrm{1}/\mathrm{2})$)./UT);
%% Eq. (11), recover flow speed along sonic path 3 before TSC

uc = [ut1.*(dT (1)./d(1));ut2.*(dT(2)./d(2));ut3.
*(dT(3)./d(3))]; %% Eq. (13)

uts1 = ut1; uts2 = ut2; uts3 = ut3;

%%Corrected 3-D wind speed

um_c=AT*uc; %% Eq. (16)

%Iteration algorithm of sonic TSC (Appendix B) for recovered data

UT_C=(um_c (1)^∧ 2+um_c ()^∧2+um_c (3)^∧2)^∧(1/2); %% Total wind magnitude

% 1st iteration

uct1 = uc (1)/(0.84+0.16.*((UT^∧2-uc (1)^∧
2)^∧(1/2))./UT);  % %flow speed 1

uct2 = uc (2)/(0.84+0.16.*((UT^∧2-uc (2)^∧
2)^∧(1/2))./UT);  % %flow speed 2

uct3 = uc (3)/(0.84+0.16.*((UT^∧2-uc (3)^∧
2)^∧(1/2))./UT);  % %flow speed 3

% 2nd iteration

for q = 2:5; %% 5 steps of iterations after 1st iteration are adequate

%TSC for flow speed 3

uct_m = [uct1(q-1);uct2(q-1);uct3(q-1)]; %% Vector of three path flow speeds

um_C = AT*uct_m; %%Vector in 3-D orthogonal system

UT_C = (um_C (1)^∧ 2+um_C (2)^∧2+um_C
(3)^∧2)^∧(1/2);
%% Total wind magnitude, again

uct3(q) = uc (3)/(0.84+0.16.*((UT_C^∧2-uct3
(q-1)^∧2)^∧(1/2))./UT_C);
%% TSC for flow speed 3

% TSC for flow speed 2

uct_mm = [uct1(q-1);uct2(q-1);uct3(q)];
%%Vector of three flow speeds, again

um_C = AT*uct_mm; %% Vector in 3-D orthogonal system, again

UT_C = (um_C (1)^∧ 2+um_C (2)^∧2+um_C
(3)^∧2)^∧(1/2); %% Recalculated the total wind magnitude

uct2(q) = uc (2)/(0.84+0.16.*((UT_C^∧2-uct2
(q-1)^∧2)^∧(1/2))./UT_C); %%% TSC for flow speed 2

%TSC for flow speed 1

uct_mm = [uct1(q-1);uct2(q);uct3(q)]; %%Vector of three flow speeds, again

um_C = AT*uct_mm; %% Vector in 3-D orthogonal system

UT_C = (um_C (1)^∧ 2+um_C (2)^∧2+um_C
(3)^∧2)^∧(1/2); %% Total wind magnitude, again

uct1(q) = u (1)/(0.84+0.16.*((UT_C^∧2-uct1
(q-1)^∧2)^∧(1/2))./UT_C); %%%TSC for flow speed 1

% Judge the steps of iterations

uct_n = [uct1(q);uct2(q);uct3(q)]; %%Vector from current iteration

ABS_C = uct_n-uct_m; %%Difference between two iterations

% Exit condition

if(abs(ABS_C(1))<  = 0.001&&abs(ABS_C(2))< 
= 0.001&&abs(ABS_C(3))<  = 0.001);

%Finalize recovered 3-D wind speed

ucm = AT*uct_n; %% Eq. (14)

ucts1 = uct1(q); ucts2 = uct2(q); ucts3 = uct3(q);

break; %% %Exit iterations

end

end

else

%Recover 3-D wind data without TSC

u = A_inversion*um; %% Acquire the flow speeds along 3 sonic paths, Eq. (10)

uc = [dT(1)./d(1).*u(1); dT(2)./d(2).*u(2);
dT(3)./d(3).*u(3)]; %%Correction

ucm = AT*uc; %%3-D orthogonal data after recovery

uts1 = uc(1); uts2 = uc(2); uts3 = uc(3);

ucts1 = ucm(1); ucts2 = ucm(2); ucts3 = ucm(3);

end

%Recover sonic temperature data

Ts = TRAW(i);

UcT =  (ucm (1)^∧2 + ucm (2)^∧2 + ucm (3)^∧2)^∧(1/2); %% Total wind

C02  =  gama_d*Rd*(Ts + 273.15); %% Eq. (32)

DELTUcT2  =  UcT^∧2 – UT^∧2;

DELTucT21  =  ucts1^∧2 – uts1^∧2; DELTucT22 = ucts2^∧2 – uts2^∧2; DELTucT23 = ucts3^∧2 – uts3^∧2;

DELTC21 = (C02 - UT^∧2 + uts1^∧ 2)*((dT(1)^∧2 –
d(1)^∧2)/d(1)^∧ 2); %% Eq. (30)

DELTC22 = (C02 - UT^∧2 + uts2^∧ 2)*((dT(2)^∧2 –
d(2)^∧2)/d(2)^∧ 2); %% Eq. (30)

DELTC23 = (C02 – UT^∧2 + uts3^∧ 2)*((dT(3)^∧2 –
d(3)^∧2)/d(3)^∧ 2); %% Eq. (30)

AAA = (DELTC21 + DELTUcT2 – DELTucT21);

BBB = (DELTC22 + DELTUcT2 – DELTucT22);

CCC = (DELTC23 + DELTUcT2 – DELTucT23);

DDD = (AAA + BBB + CCC);

EEE = 3*gama_d*Rd;

Tcs = Ts+(DDD/EEE); %% Eq. (33)

DELTC021_ad = C02*2*(1-dT(1)/d(1)); %% Eq. (38)

DELTC022_ad = C02*2*(1-dT(2)/d(2)); %% Eq. (38)

DELTC023_ad = C02*2*(1-dT(3)/d(3)); %% Eq. (38)

AAA_ad = ((dT(1)^∧2-d(1)^∧ 2)/d(1)^∧2)*(C02-
(DELTC021_ ad+((DELTC021_ad+DELTC022_ ad+
DELTC023_ad)/3))-UT^∧ 2+uts1^∧2)
+DELTUcT2-DELTucT21;

BBB_ad = ((dT(2)^∧2-d(2)^∧ 2)/d(2)^∧2)
*(C02-(DELTC022_ad+((DELTC021_ad+DELTC022_ad
+DELTC023_ad)/3))-UT^∧ 2+uts2^∧2)
+DELTUcT2-DELTucT22;

CCC_ad = ((dT(3)^∧2-d(3)^∧ 2)/d(3)^∧2)
*(C02-(DELTC023_ad+((DELTC021_ad+DELTC022_ad
+DELTC023_ad)/3))-UT^∧ 2+uts3^∧2)
+DELTUcT2-DELTucT23;

DDD_ad = (AAA_ad + BBB_ad + CCC_ad);

Tcs_ad = Ts+(DDD_ad/EEE); %% Eq. (40)

Data_recovery(i,1) =  ucm(1); %%Recovered 3-D wind speed in x-direction

Data_recovery(i,2) =  ucm(2); %% Recovered 3-D wind speed in y-direction

Data_recovery(i,3) =  ucm(3); %% Recovered 3-D wind speed in z-direction

Data_recovery(i,4) =  Tcs; %% Recovered Ts from raw Ts, Eq. (33)

Data_recovery (i,5) =  Tcs_ad; %% Recovered Ts from raw Ts, Eq. (40)

Data_recovery (i,6) =  TS_emp(i); %% Recovered T, Eq. (D1)

Data_recovery (i,7) =  TRAW(i); %% Raw Ts

End

% Output the final processing result in excel format

title = {“Recovered u_x”,“Recovered u_y”,“Recovered u_z”,“Recovered Tcs”,“Recovered Tcs_ad”,“Recovered T”,“RAW T_S”};

fname = [path “\ Data_recovery”];

xlswrite(fname,title,“sheet1”);

xlswrite(fname,Data_recovery,“sheet1”,“A2”);

waitbar(0,hwait,`Done”);

pause(2);

close(hwait);

Appendix D: Sonic temperature from air temperature, relative humidity, and atmospheric pressure

In case that air temperature (T in ^∘C), relative humidity (RH in percentage), and atmospheric pressure (P in kPa) are measured in the field, sonic temperature (T_s in ^∘C) can be calculated using the well-known equation (Kaimal and Gaynor, 1991):

$$\begin{array}{}\text{(D1)}& {\displaystyle }{T}_{\mathrm{s}}=\left(T+\mathrm{273.15}\right)\left(\mathrm{1}+\mathrm{0.32}{\displaystyle \frac{e}{P}}\right)-\mathrm{273.15},\end{array}$$

where e is air water vapor pressure (kPa) and can be computed from T, RH, and P.

Given T and P, saturated water vapor pressure (e_s in kPa) can be calculated using Buck (1981):

$$\begin{array}{}\text{(D2)}& {\displaystyle }{e}_{\mathrm{s}}=\left\{\begin{array}{l}\mathrm{0.61121}\exp \left({\displaystyle \frac{\mathrm{17.368}T}{T+\mathrm{238.88}}}\right)f_w(T,P)T\ge \mathrm{0}\\ \\ \mathrm{0.61121}\exp \left({\displaystyle \frac{\mathrm{17.966}T}{T+\mathrm{247.15}}}\right)f_w(T,P)T<\mathrm{0}\end{array}\right.,\end{array}$$

where f_w(T, P) is the enhancement factor:

$$\begin{array}{ll}{\displaystyle }{f}_w(T,P)& {\displaystyle }=\mathrm{1.00041}+P\left[\mathrm{3.48}\times {\mathrm{10}}^{-\mathrm{5}}+\mathrm{7.4}\times {\mathrm{10}}^{-\mathrm{9}}\right.\\ \text{(D3)}& {\displaystyle }& {\displaystyle }\left.{\left(T+\mathrm{30.6}-\mathrm{0.38}P\right)}^{\mathrm{2}}\right]\end{array}$$

Using the definition of air relative humidity, air water vapor pressure is given by

$$\begin{array}{}\text{(D4)}& {\displaystyle }e=e_{\mathrm{s}}{\displaystyle \frac{\mathrm{RH}}{\mathrm{100}}},\end{array}$$

Submit the measured T and P as well as the calculated e into Eq. (D1) and the sonic temperature can be calculated.