Recovery of the 3-dimensional wind and sonic temperature data from a sonic anemometer physically deformed away from manufacture geometrical settings

A sonic anemometer (sonic) reports 3-dimensional wind and sonic temperature (Ts) by measuring the time of ultrasonic signals flying along each of its three sonic paths whose geometry of lengths and angles in the sonic coordinate system was precisely determined through production calibrations and was embedded into the sonic’s firmware. If the sonic path geometry is deformed, although correctly measuring the time, the sonic continues to use its embedded geometry data for 20 internal computations, resulting in incorrect data. However, if the geometry is re-measured (i.e. recalibrated) to update sonic firmware, the sonic can resume reporting correct data. In some cases, where immediate recalibration is not possible, a deformed sonic can be used because ultrasonic signal-flying time is still correctly measured. For example, transportation of a sonic to Antarctica in 2015 resulted in a geometrically deformed sonic. Immediate deployment was critical, so the deformed sonic had been used until a replacement arrived in 2016. To recover data from this deformed sonic, equations and algorithms 25 were developed and implemented into the post-processing software to recover wind data with/without transducer shadow correction and Ts data with crosswind correction. Using two geometric datasets, production calibration and recalibration, post-processing recovered the wind and Ts data from May 2015 to January 2016. The recovery reduced the difference of 9.60 to 8.93 °C between measured and calculated Ts 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. Since such data reacquisition is time30 consuming and expensive, this data recovery approach is a cost-effective and time-saving option applicable to similar cases. The equation development can be a reference to the studies on related topics.


Introduction
The three-dimensional (3D) 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 (e.g., practically 10 to 35 50 Hz) and outputs wind speeds expressed in the 3D right-handed orthogonal anemometer coordinate system relative to its special care, back to the manufacturer of Campbell Scientific Incorporation in the US for re-measurements of its geometry to update its firmware (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 5 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.

Site, instrumentation, and data
The observation site was located 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 10 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 in particular are inhabitable to human life, but drive atmospheric dynamics in a way of interest to human beings (Valkonen et al., 2008); therefore, this region has attracted scientists to measure its surface heat balance; however, the measurements are not an easy task in financial support, technical infrastructure, and administrative management.As such, only few of studies on such measurements have been conducted in this region (e.g., Vihma et al., 15 2009;Liu et al., 2017).
Supported by National Science Foundation of China, the project: "Sea ice/snow surface energy budget of the Antarctic Prydz Bay" was initiated to measure the fluxes of CO 2 /H 2 O, heat, radiation, and atmospheric variables so that the sea ice/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 20 IRGASON (SN: 1131) for the fluxes of CO 2 /H 2 O, sensible heat, and momentum; one 4-way net radiometer (model: CNR4, Kipp & Zonen, Delft, The Netherlands) for net radiation from incoming short-wave, outgoing short-wave, incoming longwave, and outgoing long-wave radiation components; 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 25 temperature.Later 2015, a CSAT3B (Campbell Scientific Inc. UT, USA) was added for additional data of 3D wind and sonic temperature.This OPEC station is 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: Rev 04) that, in turn, was connected to and 30 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 3D wind, sonic temperature, CO 2 /H 2 O amount, atmospheric pressure, diagnosis codes for the 3D sonic anemometer and open-path infrared 35 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 every half-hour interval.Immediately after the station started to run, all measured values were checked.Unfortunately, the sonic temperature from 3D sonic anemometer was found to be incorrect because it was around 10 °C higher than the air temperature from HMP155A or 100K6A1A.Given H 2 O density about 1.00 g m -3 and air temperature about -20 °C then, 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 5 from the three sonic paths were unexpectedly various individually around -12, 5, and -7 °C 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 this difference from an IRGASON sonic anemometer was expected < 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 10 anemometer.To confirm the diagnosis, the body of IRGASON was visually examined and painting on the knuckle of side one (i.e., 1 st sonic path) among the top three claws was found off as apparently impacted (see 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 3D wind.Therefore, this IRGASON should have been shipped back to manufacturer for re-measurements of its geometry 15 to update its OS (recalibration).However, as addressed in Introduction, the 2015 data would have been missed if it were shipped back to manufacturer.To make measurements as planned, this IRGASON continued its field duty until next roundtrip of R/V Xue Long to Antarctica from China by the end of 2015 when its replacement from the manufacturer arrived at the site.

Data check and instrument diagnosis
In early 2016, it was replaced in the field and was shipped back to the manufacturer where it was re-measured for sonic 20 geometry in recalibration process on March.The re-measurements 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 need to be recovered as if measured by the same anemometer after recalibration although these data were acquired from the measurements before the recalibration.
4 Algorithm to recover the data of 3D wind and sonic temperature 25 An IRGASON sonic anemometer measures wind flows along its three non-orthogonal sonic paths (i.e. the three sonic paths non-orthogonally situated 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 to 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 flying time of ultrasonic signal upward [t ui where 30 subscript i can be 1, 2, or 3, denoting the sequential order of sonic path (see Fig. 1).This subscript denotes the same throughout] and downward (t di ) are measured by the sonic anemometer (Campbell Scientific Inc. 1998).In the case as shown in Fig. 1 for the 3 rd sonic path, or i  3, the flying time of ultrasonic signal in the path upward is given by: where, along the 3 rd sonic path, d 3 is its length precisely measured during production or recalibration process using a 35 Coordinate Measurement Machine (CMM), c 3 is the speed of sound, and u 3 is the speed of air flow (see Fig. 1); and the flying time of ultrasonic signal downward is given by: Equations ( 1) and (2) lead to: 5 Using the same procedure, u 1 and u 2 (see Fig. 1) can be derived as the same form.In reference to Eq. (3), equation for u i ; where subscript i = 1, 2, or 3; can be expressed as Similar to d 3 , d 1 and d 2 are also precisely measured using CMM.The three flow speeds of u i (i = 1, 2, or 3) measured from 10 the three non-orthogonal paths are then expressed in the 3D right-handed orthogonal instrument coordinate system of x, y, and z; where x and y are the horizontal coordinate axes and z is the vertical axis; through a transform matrix A as the 3D wind speeds (u x , u y , and u z ) commonly used in practices: where the 3D right-handed orthogonal instrument coordinate system (hereafter, sonic coordinate system.see Figs. 1 and A1) 15 is defined by its origin at the center of sonic measurement volume, the u x -u y plain parallel to the imagery plain leveled by a built-in bulb in the anemometer structure, and the u y -u z plain through the 1 st sonic path and A is a 3×3 matrix constructed using precisely measured geometry of the sonic paths in angles relative to the 3D coordinate system (see its derivations in Appendix A).Matrix A is unique for each sonic anemometer and is embedded in its firmware; therefore, the 3D wind data outputted from the anemometer are the three components of u x , u y and u z in the 3D coordinate system.20 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 IRGASON OS 5 or newer.Therefore, based on the 25 option, Eq. ( 5) alternatively can be expressed as: According to Wyngaard and Zhang (1985), the correction equation for the sonic transducer size and sonic path geometry of IRGASON sonic anemometer is given by: 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 accurately 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 where U T is the magnitude of total true wind vector, given by 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 2 nd order term under a square root in the right hand of Eq. 10 (8) although u Ti is analytically expressed in its left hand side and 2) U T is not available either because u x , u y , and u z are derived from u 1 , u 2 , and u 3 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, iteration algorithm may be a right approach to the corrections using Eq. ( 8), or to estimation of u Ti .For the 1 st iteration, u Ti in the right hand of Eq. ( 8) could be replaced with u i as its estimation.Given that U T should be 15 calculated using u x , u y , and u z from Eq. ( 6), before the 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  4) and transform matrix A in Eqs. ( 5) and ( 6) are embedded into the firmware of sonic anemometer in manufacture processes.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 25 embedded into its firmware; 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 or calibration process are erroneous.To correct the erroneous output; u x , u y , and u z need to be transformed back to t ui and t di and to 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 manufacture geometrical settings 30 before its field deployment.
For the true sonic path lengths and true sonic path angles, IRGASON (SN: 1131) was returned to the manufacturer in the way as described in Section 3. In the same way as in the manufacture process, the lengths and angles were re-measured using field uses and to correct u x , u y , u z and T s (sonic temperature, see Section.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.

i. With transducer shadow corrections
Transfer u x , u y , and u z in the 3D coordinate system to the flow speeds along the sonic paths after transducer shadow 5 corrections.
Using Eq. (B5), flow speed along the i th sonic path before transducer correction (u i ) can be expressed as 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 10 close enough to ensure u x , u y , and u z to converge at their measurement precisions (see Appendix B).Using u i and d i , the time term inside the square bracket in Eq. ( 4) can be recovered Also according to Eq. ( 4) and using d Ti , the speed of air flow along the i th sonic path can be recalculated as u ci : 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 3D wind in the 3D coordinate system ' can be expressed as: ii.Without transducer shadow corrections 20 Transfer u x , u y , and u z in the 3D coordinate system to the flow speeds along individual sonic paths Using Eqs. ( 12) and ( 13), the speed of flow along the i th 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 3D coordinate system as  1) and ( 2) also lead to: 5 Using the same procedure, c 1 and c 2 (see Figs. 1 and 5) can be derived as the same form.In reference to Eq. ( 17), equation for c i ; where subscript i = 1, 2, or 3; can be expressed as 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 10 sonic path i, is zero; c i is the true speed of sound (c Ti ).Unfortunately, crosswind rarely is zero and c i needs to be corrected to c Ti .According to Figs. 1 and 5, the true speed of sound is given by: Referencing the diagram for wind vectors in the left side of Fig. 5, this equation can be expressed as

   (20) 15
According to the definition of sonic temperature (Kaimal and Finnigan, 1994), the sonic temperature (K) along the i th sonic path (T si ) should be expressed as: where γ d (1.4003) is the ratio of dry air specific heat at constant pressure (1,004 J K -1 kg -1 ) to dry air specific heat at constant volume (717 J K -1 kg -1 ) and R d (287.04J K -1 kg -1 ) is gas constant for dry air.The sonic temperature outputted from sonic 20 anemometer (T s in °C) is the average from the three sonic paths, given by: Substituting c Ti with Eq. ( 20) and then substituting c i with Eq. (18), T s can be expressed as: Equation ( 23) indicates that, given d i , a sonic anemometer estimates sonic temperature using its measured flying 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 3D coordinate system are transformed from u i using A as explained by Eq. ( 5) without transducer corrections or from u Ti also using A as explained by Eq. ( 6) with transducer corrections.As discussed in Section 4.1.2,if a 5 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 firmware; 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 the field wind to be measured.In Section of 4.1, the procedure to recover 3D wind data was developed using re-measured sonic path lengths (d Ti ) and re-determined sonic path angles for 10 A T .The procedure to recover sonic temperature data also needs to be developed using d Ti and recovered 3D wind data in this section as follows.
Based on Eq. ( 20), the recovered speed of sound from sonic path i after crosswind corrections can be expressed as where c ci is the recovered speed of sound along sonic path i andU 2 .After replacement of c Ti 2 with 15 c cTi 2 in Eq. ( 22), the recovered sonic temperature (T cs in °C) can be written as: Now, the term of c cTi 2 needs to be derived.Subtracting Eq. ( 20) from ( 24) leads to: Using this equation to substitute c cTi 2 in Eq. ( 25), denoting U U cT T 2 2  with U cT 2 and denoting u u cTi Ti  with u cTi 2 lead to: 20 In this equation, the term of c c ci i  is still unknown.Based on Eq. ( 18), c ci 2 is given by: Accordingly, the unknown term is given by: In this equation, only unknown variable is c i 2 .Based on Eq. ( 20), this equation can be expressed as: In the right hand side of this equation, c Ti 2 is unknown only.However, the whole term in the right hand of Eq. ( 29) mathematically is a differential term in which c Ti 2 can be reasonably approximated using its neighbor value as close as possible to c Ti 2 .The average of c c c , , and can be calculated from Eq. ( 22) because T s is an output variable of sonic 5 anemometer.Without a measurement error and random error, the three c Ti 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 Ti 2 can be reasonably approximated using the average of three c Ti 2 as c T 2 , given by: where c T 2 can be computed from Eq. ( 22) as. 10  in Eq. ( 27) with Eq. ( 31) leads to 15 In the right hand side of this equation, the whole term after T s is the sonic temperature recovery term interpretable.

Application
For our case without a transducer shadow correction, Eqs. ( 15), ( 12), ( 13), and ( 16) were sequentially used to recover the 3D wind data.In a case of transducer shadow correction in option, Eqs. ( 10) to ( 16) are used.Based on the data of 3D wind from 20 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 25 Production/Calibration" whose interface is shown in Fig. 6 and Appendix C. As long as the path lengths and matrixes from production/calibration and from recalibration are input into the software as instructed by the interface (see Fig. 6), the software automatically recover the data in batches.In our station, an additional anemometer for wind was not under deployment when this studied IRGASON was used in its deformed state; therefore, no data were available to verify the recovered 3D wind data.However, the algorithms as addressed using Eq. ( 10) to Eq. ( 16) to recover the 3D wind data are solid without any estimation and the recovered 3D wind data are not necessary to be verified.

Verification
Fortunately, the data to verify sonic temperature are available in this station.Air temperature, relative humidity, and 5 atmospheric pressure were measured using research grade sensors of HMP155A and IRGASON built-in barometer and the data of these variables also stored at 10Hz (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: -20.01 ± 0.14 °C in Fig. 7a, -9.06 ± 0.13 °C in Fig. 7b, and -1.90 ± 0.22 °C in Fig. 7c.The 10 difference between measured (i.e., unrecovered) and calculated sonic temperature values of 9.60 ± 0.14 °C in Fig. 7a, 9.53 ± 0.17 °C in Fig. 7b, and 8.93 ± 0.24 °C in Fig. 7c was narrowed to 0.99 ± 0.14 °C, 0.57 ± 0.17 °C, and -0.25 ± 0.24 °C, respectively, as the difference between recovered and calculated sonic temperature values.Given the accuracy of ±0.5 °C in sonic temperature from IRGASON sonic anemometer (Personal communication with Larry Jacobsen who is the designer of sonic anemometer) and the accuracy of ±0.2 ~ 0.3 °C in air temperature below 0 °C and 1.2% in relative humidity from 15 HMP155A (Campbell Scientific Inc., 1990), 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 °C or even wider that 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 is 20 barely satisfactory in case of Fig. 7a (i.e., 0.99 ± 0.14 °C, the difference in sonic temperature between recovered and calculated) although the range of 0.99 ± 0.14 °C is not significantly different from ±0.80 °C.The bare satisfactory might be caused by the approximation of c Ti from c T that is fully valid if all c Ti are not measured by a sonic anemometer in its deformed state, but not a case in this study.
According to Eq. ( 22), it is impossible to have an individually c Ti from T s which is the sole output for sonic temperature from 25 any sonic anemometer.Now, the average of c c c , , and 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 T , , and in approximation although the exact values are impossible.The adjusted values can reflect the variability among c Ti 2 at some degree and are reasonably expected to improve the data recovery.

Adjustment 30
The measured speed of sound after crosswind correction (c Ti ) is independent of wind speed (Liu et al. 2001).Given air density and atmospheric pressure (Barrett and Suomi, 1949), without wind, c Ti 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 Given the flying time is correctly measured by a sonic anemometer (i.e., 0 i t ) even in its geometrical deformation, this equation becomes Using this equation to replace c Ti 2 in Eq. ( 30) and the resultant equation with this replacement then is used to c c ci i 2 2  in Eq.
(27) as In the right hand side of equation, the whole term after T s is the adjusted sonic temperature recovery term.The data ever recovered using Eq. ( 33) also were recovered again using Eq.(40).Apparently, this equation did a better job 25 than Eq. ( 33).The difference in sonic temperature between recovered and calculated was reduced to 0.81 ± 0.14 °C, 0.38 ± 0.17 °C, and -0.45 ± 0.24 °C, respectively, as shown from panels a to c in Fig. 7.These values for the difference fall into the  40) is believed to do a better job than Eq. ( 33), although, that is satisfactory.Eventually, Eq. ( 40) was used for data recovery and was incorporated into the software as shown in Fig. 6 and Appendix D.

Verification of 3D wind recovery
Although not explicitly verified, the recovered 3D wind data were implicitly verified through the verification of recovered 5 sonic temperature data because the recovery of sonic temperature data must rely on recovered 3D wind data [see Eqs. ( 33) and ( 40)] as the cross wind correction for sonic temperature needs 3D wind data (Liu et al., 2001).If 3D wind had not been well recovered, sonic temperature data could not have been recovered satisfactorily.The satisfactory recovery of sonic temperature data in this study implicitly verified the satisfactory recovery of 3D wind data.

Comparability of recovered to calculated sonic temperature 10
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 IRGASON built-in barometer (see Appendix D); Therefore, the former reflected the fluctuations in sonic temperature at high frequency, the latter reflects such fluctuations not so high as the former.As such, a pair of recovered and calculated sonic temperature values from simultaneous measurements (i.e., the same record in a time series data file) were 15 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.
For better comparison, the difference was calculated from the recovered and calculated sonic temperature values temporally aligned in consideration of lag.The calculated sonic temperature from the air temperature, relative humidity, and atmospheric pressure; which were measured using slow response sensors; was believed to be lagged behind recovered one in 20 response to the fluctuations in sonic temperature.The lag time about 10 seconds was empirically found at the maximization of cross correlation of a time series of recovered sonic temperature to different time-lagged (i.e., time-shifted) series of calculated sonic temperature (Ibrom et al., 2007).The difference between recovered and calculated sonic temperature values was calculated using recovered sonic temperature without a time lag and calculated sonic temperature with a time lag of 10 seconds.The use of lag time might be unnecessary, but might make the comparison as reasonable as possible.25 8.3 Recovered higher than calculated sonic temperature at lower temperature See Fig. 7. Compared to calculated sonic temperature, the recovered sonic temperature from Eq. ( 40) is 0.81 ± 0.14 °C higher at -20.01 °C (Fig. 7a) and 0.38 ± 0.17 °C higher at -9.06 °C (Fig. 7b), however, at -1.90 °C, even 0.45 ± 0.24 °C 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 (Campbell 30 Scientific Inc., 1990).
The sonic path lengths and geometry of sonic anemometer were measured at the manufacture environment of air temperature around 20 °C (i.e., manufacture temperature) and embedded into its firmware 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 manufacture temperature while the length values of sonic paths are 35 unchanged inside firmware.As a result, the sonic anemometer could under-or over-estimate the speed of sound, thus sonic temperature.The under-or over-estimation 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 the case of this study, the working air temperature for the sonic anemometer was as low as negative to -20 °C, within which the sonic paths become shorter at some degree so that its measurement performance possibly was 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 centigrade of tens below manufacture temperature and thus 5 sonic temperature is anticipated as shown in Fig. 7a to Fig. 7c.
Although, at different air temperature, the performance of the temperature and relative humidity probe and IRGASON builtin barometer, whose measurements are used to calculate the sonic temperature (see Appendix D), is more stable than a sonic anemometer although their accuracies are the best at 20 °C, too, and become lower with temperature away from 20 °C (Campbell Scientific Inc., 1990).For example, HMP155A has an accuracy in air temperature to be ±0.1 °C at 20 °C and to 10 be ±0.25 °C at -20 °C as well as an accuracy in relative humidity (RH) to be ±(1.0+0.008RH) in % at 20 °C and to be ±(1.2+0.012RH) in % at -20 °C.The greater disagreement between recovered and calculated sonic temperature values at lower temperature in Fig. 7a may also be contributed by the fact that the lower the air temperature, the lower the accuracies of HMP155A and IRGASON built-in barometer.

Radiation on calculated sonic temperature 15
See Fig. 7c.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 in part by higher incoming solar radiation of 750 W m -2 in short-wave on the radiation shield of HMP155A.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 IRGASON built-in barometer.The 20 HMP155A housed inside a radiation shield (see Fig. 2) was subject to contamination from solar radiation.Even 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., 25 2009).In case of Fig. 7c at 750 W m -2 of incoming short-wave radiation, a degree higher of air 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 (see Eq. D1 in Appendix D)

Applicability of equations and algorithms in this study
Any sonic anemometer is built as slender (e.g., < 1.00 cm in each diameter of CSAT3 six claws to hold individual sonic 30 transducers) and 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 in transportation (e.g., the case in this study), installation, or other handlings.As shown in this study, a slight geometrical deformation, even changes of millimeters or less of sonic path length (see Table A1 in Appendix A) could cause significant errors in 3D wind and, especially, in sonic 35 temperature.According to our recalibration experience with 3D sonic anemometers at Campbell Scientific Incorporation, these cases as addressed in this study have been not unusual, but the equations and algorithms to recover the data measured by a deformed 3D sonic anemometer were not available.Since requisition of these datasets are expensive, thus their recovery would be the 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 IRGASON sonic anemometer.The physics is the same as those for other models of Campbell Scientific 3D sonic anemometers such as CSAT3, CSAT3A, and CSAT3B that are popularly used in the world (Horst et al., 2015).The sonic path geometry of IRGASON sonic anemometer, however, is different from other models in the assigned azimuth angle of 1 st sonic path in the 3D coordinate system.This angle was assigned as 90° in IRGASON sonic anemometer, but as 0° in 5other models (Campbell Scientific Inc., 1998;2010b;2015).Even so, given the sonic path lengths and transfer matrixes of sonic anemometer that were measured and determined in manufacture or calibration process [d i in Eq. ( 12) and A in Eq. ( 15)] and in recalibration process after the use in its 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 3D sonic anemometers (see Fig. 6).The derivation procedures and even equations based on the measurement working physics are 10 applicable as a reference to the development of the equations and algorithms to recover the data measured using other brands of 3D sonic anemometers that incurred deformations or to studies on similar topics.

Conclusion remarks
An IRGASON 3D sonic anemometer (SN: 1131) which was geometrically deformed by impact during transportation to Antarctica from China early 2015.To fulfill the field measurement plans for the year, it had to be deployed there in the 15 Zhongshan Station until early 2016 when it was replaced in the field with another IRGASON provided by the manufacturer and was returned to its manufacturer, Campbell Scientific Incorporation, for recalibration through the re-measurements of its sonic path geometry (i.e., lengths and angles), re-determination of transfer matrix, and update of operating system.To recover the 3D 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 20 IRGASON/CSAT3/A/B Used in Geometrical Deformation after Production/Calibration" as shown in Fig. 6.Given two sets of sonic path lengths and two transfer matrixes of sonic anemometer that were measured and determined in manufacture/calibration process and also in recalibration process after the use in its deformed state, the data measured by the IRGASON 3D sonic anemometer even in its deformed state were recovered as if measured by the same anemometer recalibrated immediately after its deformation.25 Inside a Campbell Scientific sonic anemometer, the transducer shadow correction for 3D wind (Kaimal and Finnigan, 1994) is a programmable option to a user; however, the crosswind correction for sonic temperature (Liu et al., 2001) is internally applied as default by its firmware.In a case of transducer shadow correction in option, the 3D 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 3D wind, the sonic temperature data are recovered using Eqs.( 9), ( 32), (38), and (40); therefore, the satisfactory recovery for 30 both 3D wind data and sonic temperature can be reflected eventually by the satisfactory of sonic temperature data recovery.
The software based on the equations and algorithms from this study can recover the 3D wind data with/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 (see Appendix D).As shown in Fig. 7, the recovered data of sonic temperature using Eq. ( 33) and Eq. ( 40) were compared to the calculated sonic 35 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 °C between unrecovered and calculated sonic temperature (i.e., unrecovered minus calculated) was narrowed by Eq. ( 40) to 0.81 to -0.45 °C (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) 3D sonic anemometer in its deformed state from May 2015 to January 2016.The eight-month data were recovered using three days of one engineer's time.Further using EddyPro 6.2.0 (Li-Cor Inc., 2016), the recovered data were further processed for the fluxes of CO 2 /H 2 O, sensible heat, and momentum.The data quality (Foken et al., 2012) mostly range 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 -1 .5 The use of a deformed 3D sonic anemometer is a practical case.If the data from such a use cannot be recovered, the requisition of these data are 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 3D sonic anemometers such as CSAT3, CSAT3A, and CSAT3B that are popularly 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 10 and algorithms to manage the data measured using other brands of 3D 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 (3D) orthogonal coordinate 15 system of instrument or natural wind, but a sonic anemometer measures flow velocities along its three non-orthogonal sonic paths (i.e.situated non-orthogonally each other, see Figs. 1 and A1); therefore, for applications, the flow velocities along the three sonic paths need to be transformed into a 3D right-handed orthogonal coordinate system in reference to the geometry of sonic anemometer as shown in Fig. A1 (i.e., the 3D orthogonal instrument 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 20 are the three coordinate axes in the 3D orthogonal instrument coordinate system.This system is defined with the x-y plain parallel to the anemometer bulb-leveled instrument plain, with the 1 st sonic path on the y-z plain, and with origin in the center of measurement volume.A flow speed along the ith (i = 1, 2, or 3) sonic path is a combination of component velocities along the path from u x , u y , and u z ; given by: where θ i and φ i are the zenith and azimuth angles of the ith sonic path in the 3D orthogonal instrument coordinate system.In this system (see Fig. A1), given the 1 st sonic path has an azimuth angle of φ 1 equal to 90° as fixed on the x-y plain, Eq. (A1) can be expressed in a matrix form of sin cos sin cos sin sin cos sin cos sin sin cos where A is a matrix expressing the flow speeds along the three non-orthogonal sonic paths in the 3D orthogonal instrument 30 coordinate system.Nominally for the sonic paths of IRGASON, θ 1 , θ 2 , and θ 3 are all 30° and φ 2 and φ 3 are 330° and 210° (see Fig. A1).Given φ 1 = 90°, these angles are calculated using measured data from Coordinate Measurement Machine and, along with the sonic path lengths, are listed in Table A1 for IRGASON Serial Number of 1131 before and after its geometrical deformation.T were used for our data recovery and A was also used in the firmware inside the IRGASON sonic 15 anemometer.

Appendix B Iteration algorithm for sonic transducer shadow corrections
Given transform matrix A, using Eq. ( 5), the measured wind vector   .Subsequently, U T is calculated using Eq. ( 9).Replace u Ti with u i under the square root in the right hand of Eq. ( 8), an approximate equation for 20 the 1 st iteration is given: where subscript i is 1, 2 or 3 and subscript _1 of u Ti indicates that it is calculated from the 1 st iteration.

st iteration
Equation ( B1) is used for sonic transducer shadow corrections in the first iteration.
Using Eq. ( 9), U T is recalculated.Replace u i with u Ti_1 under the square root in the right hand of Eq. ( B1), an approximate equation for the 2 nd iteration is given: .
Given T and P, saturated water vapor pressure (e s in kPa) can be calculated using Buck ( 1981): Submit the measured T and P as well as the calculated e into Eq.(D1), the sonic temperature can be calculated.
relative error of whole term in the right hand side of Eq. (31) would be < 4% even if the variability in sonic temperature due to the difference among c Ti 2 values reaches 10 °C at air temperature of -30 °C without wind (i.e.,UT  0 and u Ti  0 ), which would be the worst case.Substituting the term of c Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2018-92Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 22 June 2018 c Author(s) 2018.CC BY 4.0 License.
34) 35 Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2018-92Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 22 June 2018 c Author(s) 2018.CC BY 4.0 License.12 In Eq. (33), c T 2 is the average of three squared c Ti [see Eqs.(22) and (32)], but an individual c Ti is unknown; therefore, for recovery improvement, it has to be estimated from c T 2 through a reasonable adjustment.The difference in magnitude between added or subtracted with the same constant to make the average of three adjusted c Ti 2 values as c T 2 , but the variability among c Ti 2 values is kept the same.This constant must be the mean of three c Ti Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2018-92Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 22 June 2018 c Author(s) 2018.CC BY 4.0 License.range of ±0.80 °C in statistical sense.Equation ( Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2018-92Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 22 June 2018 c Author(s) 2018.CC BY 4.0 License.
subscript T indicates "True" because, after IRGASON deformation, it should be used in the field although it was not used.The inversion of this matrix is given as Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2018-92Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 22 June 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 1
Figure 1 Diagram of IRGASON for the three sonic measurement paths (red dash lines) along which ultrasonic signals fly and the three dimensional (3D) right-handed orthogonal instrument coordinate system (blue lines) in which 3D wind is expressed (i.e.u 1 , u 2 , and u 3 are the flow speeds along the 1 st , 2 nd , and 3 rd sonic paths, respectively.These three flow speeds are expressed as u x , u y , 5

Figure 4 :
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 i th sonic path; u x , u y , and u z are the wind speeds expressed in the three-dimensional orthogonal instrument coordinate system; and α i is the angle between sonic path i and the total flow vector (U T ) equal to u u i i 2 2   or u u u x y z 2 2 2   ].See Wyngaard and Zhang (1985) and Kaimal and

Figure 5 :
Figure 5: Crosswind on 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 u u i i 2 2   or u u u x y z 2 2 2 

Figure 7 : 5 [
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 IRGASON built-in barometer { T s measured by the IRGASON sonic anemometer in geometrical deformation (raw T s ), T s recovered from raw T s using equation (33), T s recovered also from raw T s using equation (40) 5

Figure A1 :
Figure A1: IRGASON sonic path angle geometry in the three-dimensional right-handed instrument 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, a red arrow below the IRGASON is the projection of the corresponding sonic path vector on the x-y plain, i.e. instrument bubble-leveled plain.As indicated by their subscript of 1, 2, or 3 for the 1 st , 2 nd , or 3 rd sonic path, θ 1 , θ 2 , and θ 3 are 5

Table A1：The lengths, zenith angles, and azimuth angles of sonic paths in IRGASON (Serial Number: 1131) instrument coordinate system before and after its geometrical deformation (measured using Coordinate Measurement Machine in September 09, 2014 before the deformation and in March 06, 2016 after use in deformation)
Using the data in this table, matrix A in Eq. (A2) and its inversion A -1 for this IRGASON before its geometric deformation 5 (i.e., as used in IRGASON firmware although not valid in the field after deformation) are given