## Methodology
\begin{bmatrix}u\\ v\end{bmatrix} _{p}\] (2) \[= \begin{bmatrix}\sum_{k=1}^{K}(\alpha_{k}w_{k}v_{LOS_{k}})\\ \sum_{k=1}^{K}(\beta_{k}w_{k}v_{LOS_{k}})\end{bmatrix} _{p},\]in which \(w_{k}\) is a weighting parameter, expressed as: \[w_{k}=\frac{R^{2}-r_{k}^{2}}{R^{2}+r_{k}^{2}}\cdot\frac{K}{K_{l_{k}}L},\] (3) where \(r_{k}\) is the distance between each measurement \(k\) to the grid point \(p\), \(K_{l_{k}}\) is the number of measurements contributed by the lidar \(l\) which corresponds to measurement \(k\) and \(L\) is the total count of contributing lidars. Different weighting functions have been used by other researchers for similar purposes [25; 26], including the range gate weighting function of pulsed lidars [27]. The Cressman weighting procedure [25] was validated by comparison with sonic anemometer measurements by [PERSON] et al. [28]. In Equation (3), the first fraction resembles this Cressman weight, which is a parabolic function that assigns measurements a weight based on their distance \(r_{k}\) to a grid point \(p\) inside the circle with radius of influence \(R\). The second fraction is the lidar weight assigned to each single measurement \(k\), which equalises the contributions from the lidars to the grid point. An additional weight considering the accuracy of each line-of-sight measurement, e.g., based on Carrier-to-Noise-Ratio (CNR), could also be included, but this is not addressed in this paper. 4. For a grid of size \(m\)-by-\(n\), there will be a number of \(mn\) grid points \(p\). Instead of solving Equation (2) for each grid point \(p\), a sparse \(2 mn\)-by-\(2 mn\) matrix \(\mathbf{A}\) is defined. Each of the four elements inside the matrix in Equation (2) will form a diagonal with length \(mn\) on the corresponding quadrant of the newly defined sparse matrix \(\mathbf{A}\). Now, the larger set of linear equations is defined: \[\mathbf{A}\,\mathbf{x}=\mathbf{b},\] (4) with: \[\mathbf{x}= \begin{bmatrix}u_{1}\\ \vdots\\ \vdots\\ u_{mn}\\ v_{1}\\ \vdots\\ v_{mm}\end{bmatrix} ,\quad\mathbf{b}= \begin{bmatrix}\sum_{k=1}^{K_{1}}(\alpha_ {k}w_{k}v_{LOS_{k}})_{1}\\ \vdots\\ \vdots\\ \sum_{k=1}^{K_{mn}}(\alpha_{k}w_{k}v_{LOS_{k}})_{mn}\\ \sum_{k=1}^{K_{1}}(\beta_{k}w_{k}v_{LOS_{k}})_{1}\\ \vdots\\ \sum_{k=1}^{K_{mn}}(\beta_{k}w_{k}v_{LOS_{k}})_{mm}\end{bmatrix} .\] (5)
5. The continuity adjustment is included in the linear system in the form of additional rows, in order to establish an _over-determined_ linear system. This provides a more robust system in case \(\mathbf{A}\) is ill-defined, which could occur e.g., when there are local scarcities in the measurements or when a part of the line-of-sight measurements is linearly dependent. The two-dimensional continuity equation for incompressible flow can be written in differential form as: \[\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0.\] (6) For inclusion in the system, a numerical approximation of the continuity equation is necessary. Where applicable, the central difference scheme is used for expressing the partial differentials. At the grid boundaries, either upward or downward difference schemes [29] are needed due to a lack of neighboring points. By adding the continuity adjustment to \(\mathbf{A}\), a new system with the \(3 mn\)-by-\(2 mn\) matrix \(\mathbf{A_{C}}\) is created. There are still \(2 mn\) unknown variables, leaving the vector \(\mathbf{x}\) unaffected. The vector \(\mathbf{b_{C}}\) is generated by extending \(\mathbf{b}\) with \(mn\) zeros. The final linear system,  
\[\mathbf{A_{C}}\,\mathbf{x}=\mathbf{b_{C}}, \tag{7}\]  
is an over-determined linear system, which is solved by minimising its squared residuals.