2.1. State space model

In this section, the 4D-Var problem in the context of turbulent flow reconstruction is formulated. The problem aims to reconstruct the turbulent flow state using a time series of lidar measurements and an LES model. As shown by Stuart (Reference Stuart2010), the continuous 4D-Var problem needs to be formally derived in a functional space, requiring us to deal with probability measures with infinite dimensions. In order to avoid this complexity, Lorenc (Reference Lorenc1986) suggested to represent the model in a finite discrete basis before deriving the 4D-Var problem. For instance, Bauweraerts & Meyers (Reference Bauweraerts and Meyers2020) formulated their 4D-Var problem in a truncated KL basis, which led to a mathematically simpler problem formulation. In the current study, we exploit the spatial discretization of the governing equation (see Appendix A) to derive the strong 4D-Var problem directly in a discrete basis. This is further elaborated below.

Consider a velocity vector $\boldsymbol {U}_0=[U_0^1,U_0^2,U_0^3]^T\in \mathbb {R}^N$ defined at time $t=0$ on a domain of size $L_1 \times L_2 \times H$, with $N_1 \times N_2 \times N_3$ uniformly distributed grid points. Our current discretization uses a staggered grid for the vertical velocity component (see Appendix A for details), so that $N=N_1N_2(3N_3-1)$. Given the initial condition $\boldsymbol {U}_0$, the LES model can be used to progress the solution in time. We write in short

(2.1)\begin{equation} {{\boldsymbol{U}}}(t)=\mathcal{M}_t(\boldsymbol{U}_0),\end{equation}

where $\mathcal {M}_t$ is the solution operator representing the time integration of the LES model described in Appendix A and ${{\boldsymbol {U}}}(t)$ is the (modelled) state at time $t$. Furthermore, the initial condition $\boldsymbol {U}_0$ needs to be solenoidal. This leads to a constraint in the 4D-Var problem, i.e. the discretized divergence of the initial velocity $\boldsymbol {U}_0$ must be zero (see Appendix C). To avoid a constrained optimization problem, we simply eliminate the continuity constraint. To this end, we construct a projection and reconstruction operator such that $\boldsymbol {\varTheta }_0=\mathcal {P}(\boldsymbol {U}_0)\in \mathbb {R}^M$ with $M=N_1N_2(2N_3-1)$, $\boldsymbol {U}_0=\mathcal {R}({\boldsymbol {\varTheta }_0})$ with $\boldsymbol {\varTheta }_0=\mathcal {P}(\mathcal {R}(\boldsymbol {\varTheta }_0))$ and $\text {div}(\mathcal {R}({\boldsymbol {\varTheta }_0}))=0$. More details are provided in Appendix C.

2.2. Methodology

In this section, the strong 4D-Var problem is introduced. Consider a series of $N_s+1$ observations ${\boldsymbol{\mathsf{Y}}} \triangleq [\boldsymbol {y}_0, \ldots, \boldsymbol {y}_{N_s}]$ sampled every $T_s$, with $\boldsymbol {y}_k \in \mathbb {R}^{N_{m}}$ a vector of $N_{m}$ measurements at time $t_k$, where $t_k=t_0+kT_s$. The latter vector can be modelled as

(2.2)\begin{equation} \boldsymbol{y}_k=\boldsymbol{h}_k(\mathcal{M}_k(\boldsymbol{U}_0)+ \boldsymbol{\epsilon}_k)+\boldsymbol{v}_k,\end{equation}

where $\boldsymbol {h}$ is the observation function introduced in § 5 for the lidar measurements and ${\boldsymbol {v}}$ and ${\boldsymbol {\epsilon }}$ are the measurement and model errors, respectively. Using the strong assumption that $\boldsymbol {h}_k(\mathcal {M}_k(\boldsymbol {U}_0)+\boldsymbol {\epsilon }_k)=\boldsymbol {h}_k(\mathcal {M}_k(\boldsymbol {U}_0))+\boldsymbol {h}_k(\boldsymbol {\epsilon }_k)$, which is valid for linear observation functions, and assuming that $\boldsymbol {\xi }_k=\boldsymbol {h}_k(\boldsymbol {\epsilon }_k)+\boldsymbol {v}_k$ is independent and normally distributed with variance $\gamma ^2\boldsymbol {I}$, we can express Bayes’ rule as $p(\boldsymbol {\varTheta }_0\mid {\boldsymbol{\mathsf{Y}}})\sim p({\boldsymbol{\mathsf{Y}}}\mid \boldsymbol {\varTheta }_0) p(\boldsymbol {\varTheta }_0),$ with

(2.3)\begin{equation} p({\boldsymbol{\mathsf{Y}}}\mid\boldsymbol{\varTheta}_0) \propto \prod_{k=1}^{N_s} \exp (-{\|\boldsymbol{y}_k-\boldsymbol{h}_k(\mathcal{M}_k (\mathcal{R}({\boldsymbol{\varTheta}_0})))\|^2}/{2 \gamma^2}) \end{equation}

and $p(\boldsymbol {\varTheta }_0)\propto \exp (-{{\boldsymbol {\varTheta }}}_{0}^{T} \boldsymbol {B}^{-1}{{\boldsymbol {\varTheta }}}_0/2)$ with ${\boldsymbol {B}}={\mathcal {P}} {\boldsymbol {R}} {\mathcal {P}}^T$, where $\boldsymbol {R}$ is the two-point correlation tensor defined in § 3. In order to obtain the best estimate of the state $\boldsymbol {\varTheta }_0$, the MAP estimate $\arg \max _{\boldsymbol {\varTheta }_0} p({\boldsymbol{\mathsf{Y}}}\mid \boldsymbol {\varTheta }_0) p(\boldsymbol {\varTheta }_0)$ is often used. This can be obtained by taking the logarithm as $-\log [p(\boldsymbol {\varTheta }_0\mid {\boldsymbol{\mathsf{Y}}})]$, leading to the following unconstrained optimization problem over the initial state $\boldsymbol {\varTheta }_{0}$:

(2.4)\begin{equation} \underset{{{\boldsymbol{\varTheta}}}_{0}}{\operatorname{minimise}} \quad \mathscr{J}({{\boldsymbol{\varTheta}}}_0)=\underbrace{\frac{1}{2} {{\boldsymbol{\varTheta}}}_{0}^{T} \boldsymbol{B}^{{-}1} {{\boldsymbol{\varTheta}}}_0}_{\mathscr{J}_B} + \underbrace{\frac{1}{2 \gamma^2} \sum_{n=1}^{N_s}\| \boldsymbol{y}_n-\boldsymbol{h}_n(\mathcal{M}_t({\mathcal{R}}( {{\boldsymbol{\varTheta}}}_0)) \|^2}_{\mathscr{J}_L}.\end{equation}

In the current study, our focus is on the initial background tensor $\boldsymbol {B}$ and the regularization term $\mathscr {J}_B$. The latter term plays a crucial role in ensuring a well-posed assimilation problem and promoting fast convergence. Furthermore, it holds particular significance in regularizing the solution in areas where no measurements are available. At one extreme, the most straightforward approach for the covariance tensor is to use the identity matrix, $\boldsymbol {R}=\boldsymbol {I}$ (or $\boldsymbol {B}={\mathcal {P}}{\mathcal {P}}^T$), which is equivalent to the standard ridge regression. However, this approach leads to slow convergence and relatively high reconstruction errors, as will be further discussed. At the other extreme, an ideal approach would involve numerically constructing the true background tensor, which was demonstrated by Bauweraerts & Meyers (Reference Bauweraerts and Meyers2020) to be prohibitively expensive.

In an attempt to achieve a trade-off between accuracy and simplicity, we study the feasibility of replacing the true background tensor with simple yet realistic analytical approximations. To this end, we use two different models for the two-point covariance tensor $\boldsymbol {R}$, which are discussed in detail in the following section.

3.1. Overview

For a given mean velocity field $\bar {\boldsymbol {u}}(\boldsymbol {x})$, the fluctuations ${\boldsymbol {u}}^{\prime }(\boldsymbol {x})= \boldsymbol {u}(\boldsymbol {x})-\bar {\boldsymbol {u}}(\boldsymbol {x})$ can be defined. The two-point velocity-correlation tensor is defined as

(3.1)\begin{equation} R_{i j}(\boldsymbol{x},\boldsymbol{\breve{x}})=\langle u_i^{\prime}(\boldsymbol{x}) u_j^{\prime}(\boldsymbol{\breve{x}})\rangle,\end{equation}

where $\langle {\cdot } \rangle$ indicates ensemble averaging. For homogeneous turbulence, the correlation tensor depends only on the separation $\boldsymbol {r}=(\breve {x}_1-x_1, \breve {x}_2-x_2, \breve {x}_3-x_3)$ between the two points. Hence, the isotropic–homogeneous spectral density tensor $\hat {\boldsymbol {\varPhi }}^{iso}(\boldsymbol {k})$, with wavenumber vector $\boldsymbol {k}=(k_1,k_2,k_3)$, can be obtained by applying the Fourier transform

(3.2)\begin{equation} \hat{\boldsymbol{\varPhi}}^{iso}(\boldsymbol{k}) = \frac{1}{8 {\rm \pi}^3} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} R_{i j}(\boldsymbol{r}) \\ \exp(-\textrm{i}(k_1 r_1+k_2 r_2+k_3 r_3))\, \mathrm{d} r_1\, \mathrm{d} r_2 \,\mathrm{d} r_3.\end{equation}

For incompressible flows, the latter tensor can be expressed in terms of the energy spectrum $E(k)$ (Pope Reference Pope2000) as

(3.3)\begin{equation} \hat{\varPhi}_{ij}^{iso}(\boldsymbol{k} )=\frac{E(k)}{4 {\rm \pi}k^4}(\delta_{i j} k^2-k_i k_j), \end{equation}

with $k=(k_1^2+k_2^2+k_3^2)^{1/2}$. An expression for the energy spectrum $E(k)$ can be analytically derived (e.g. Pope Reference Pope2000) based on the asymptotic behaviour of the turbulent energy for large and small scales. In this paper, we investigate two variants of the energy spectrum $E(k)$ based on the asymptotic behaviour as $k \rightarrow 0$. The spectrum can be expressed as

(3.4)\begin{equation} E(k)=a \sigma_{iso}^2 \ell \frac{(k\ell)^p}{(1+(k\ell)^2)^{{5/6}+{p/2}}},\end{equation}

where $\ell$ is a length scale and $\sigma _{iso}^2$ is the isotropic variance. The power $p$ determines the slope of $E(k)$ as $k \rightarrow 0$. When $p=4$, the energy spectrum exhibits a behaviour of $E(k) \sim k^4$ for large scales, commonly known as Batchelor turbulence (Davidson Reference Davidson2010). Conversely, $p=2$ leads to a behaviour of $E(k) \sim k^2$, which is also known as Saffman turbulence (Saffman Reference Saffman1967; Pope Reference Pope2000). The scaling factor $a$ is chosen such that the integration of $E(k)$ over all wavenumbers yields the total turbulent energy $E_{tot}$ for both Batchelor and Saffman cases. This results in values of 1.45 and 1.1886 for the Batchelor and Saffman cases, respectively.

Although assuming homogeneous isotropic turbulence leads to great mathematical simplifications, wall-bounded flows in reality are known to behave otherwise. Nevertheless, in a neutral pressure-driven ABL over smooth terrain, turbulence is close to homogeneous in the horizontal plane (Holmes, Lumley & Berkooz Reference Holmes, Lumley and Berkooz1996), and therefore homogeneity may be assumed in the $x_1$ and $x_2$ directions but not in the $x_3$ direction. Consequently, the two-point correlation functions become dependent on the absolute vertical location of each point rather than the distance between them, and the Fourier transform over the vertical distance in (3.2) may not be applied. Instead, the two-point correlations are expressed in terms of the vertically inhomogeneous spectral tensor $\hat {\boldsymbol {R}}$, which can be constructed from the function

(3.5)\begin{align} \hat{R}_{ij}(k_1,k_2;x_3, \breve{x}_3)=\frac{1}{4 {\rm \pi}^2} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} R_{i j}(r_1, r_2;x_3, \breve{x_3}) \exp (-\textrm{i}(k_1 r_1+k_2 r_2))\, \mathrm{d} r_1\, \mathrm{d} r_2.\end{align}

Owing to our pseudo-spectral numerical discretization (see Appendix A), we focus on obtaining the correlation function directly as $\hat {R}_{ij}(k_1,k_2;x_3, \breve {x}_3)$ instead of $R_{i j}(\boldsymbol {x},\boldsymbol {\breve {x}})$ in the remainder of this section.

3.2. The HGW model

Hunt & Graham (Reference Hunt and Graham1978) studied the effect of introducing a moving rigid surface suddenly at time $t=0$ on the isotropic homogeneous tensor ${\boldsymbol {\varPhi }}^{iso} (\boldsymbol {k})$. To model the blockage by the wall, the velocity field $\boldsymbol {u}(\boldsymbol {x},t)$ was written as a sum of the homogeneous flow field $\boldsymbol {u}^{H}(\boldsymbol {x},t)$ from the homogeneous tensor and a blockage field $\boldsymbol {u}^{B}(\boldsymbol {x},t)$ obeying $u_3^{H}(x_1,x_2,0,t)+u_3^{B}(x_1,x_2,0,t)=0$. When the shear-free case is considered (i.e. zero mean shear), $\boldsymbol {u}^{B}$ is simply an irrotational velocity field that does not affect the initial vorticity field (Lee & Hunt Reference Lee and Hunt1989). For further details, we refer the reader to the original paper. Wilson (Reference Wilson1997) used these ideas to derive an analytical model for the vertically inhomogeneous spectral tensor $\hat {\boldsymbol {R}}(k_1,k_2;x_3, \breve {x}_3)$ in the atmospheric CBL. The model was named the HGW model and is given as

(3.6)\begin{align} \hat{R}_{ij}^{HGW} (k_1,k_2;x_3, \breve{x}_3) &= \hat{R}_{ij}^{iso} (k_1, k_2 ;\breve{x}_3-x_3) \nonumber\\ &\quad +{\rm e}^{-\kappa \breve{x_3}} m_j(k_1, k_2) \hat{R}_{3 i}^{{iso}*}(k_1, k_2;x_3) \nonumber\\ &\quad +{\rm e}^{-\kappa x_3} m_i^*(k_1, k_2) \hat{R}_{3 j}^{iso}(k_1, k_2;\breve{x}_3) \nonumber\\ &\quad +{\rm e}^{-\kappa(x_3+\breve{x_3})} m_i^*(k_1, k_2) m_j(k_1, k_2) \hat{R}_{33}^{iso}(k_1, k_2;0), \end{align}

where $\kappa =(k_1^2+k_2^2)^{1/2}$ and $(m_1, m_2, m_3)=(\textrm {i} k_1 / \kappa, \textrm {i} k_2 / \kappa,-1)$, with $\hat {R}_{ij}^{iso} (k_1, k_2, r)$ obtained by computing the inverse Fourier transform (IFT) of the three-dimensional spectral tensor $\hat {\varPhi }^{iso}_{ij} (\boldsymbol {k})$ in the direction of $k_3$, which can be analytically computed for Batchelor turbulence (Wilson Reference Wilson1997) (see Appendix D for more details). Although the model in (3.6) accounts for the vertical inhomogeneity, it does not capture the anisotropic behaviour imposed by the mean flow. The implications of this are further discussed in § 6.

3.3. The Mann model

Many turbulent tensor models in the literature are inherently isotropic. However, anisotropy effects are often included by introducing stretching factors in different directions, leading to many extra coefficients to be tuned (in addition to the isotropic tensor parameters) (Kristensen et al. Reference Kristensen, Lenschow, Kirkegaard and Courtney1989). Alternatively, Mann (Reference Mann1994) introduced a turbulence model based on rapid distortion theory, which accounts for anisotropic effects via a single additional parameter $\varGamma$. In this model, the spectral tensor is assumed to be initially isotropic at time $t=0$, as given in (3.3). Afterwards, the isotropic tensor is distorted by a constant shear profile for a certain period of time that is equal to the eddy lifetime. The resulting model is summarized as

(3.7a)$$\begin{gather} \hat{\varPhi}_{11}^{Mann}(\boldsymbol{k}) = {E (k_0 )} (k_0^2-k_1^2-2 k_1 k_{30} \zeta_1+ (k_1^2+k_2^2 ) \zeta_1^2 )/({4 {\rm \pi}k_0^4}), \end{gather}$$

(3.7b)$$\begin{gather}\hat{\varPhi}_{22}^{Mann}(\boldsymbol{k}) = {E (k_0 )} (k_0^2-k_2^2-2 k_2 k_{30} \zeta_2+ (k_1^2+k_2^2 ) \zeta_2^2 )/({4 {\rm \pi}k_0^4}), \end{gather}$$

(3.7c)$$\begin{gather}\hat{\varPhi}_{33}^{Mann}(\boldsymbol{k}) = {E (k_0 )} (k_1^2+k_2^2 )/({4 {\rm \pi}k^4}), \end{gather}$$

(3.7d)$$\begin{gather}\hat\varPhi_{12}^{Mann}(\boldsymbol{k}) = {E (k_0 )} ({-}k_1 k_2-k_1 k_{30} \zeta_2-k_2 k_{30} \zeta_1+ (k_1^2+k_2^2 ) \zeta_1 \zeta_2 )/({4 {\rm \pi}k_0^4}), \end{gather}$$

(3.7e)$$\begin{gather}\hat\varPhi_{13}^{Mann}(\boldsymbol{k}) = {E (k_0 )} ({-}k_1 k_{30}+ (k_1^2+k_2^2 ) \zeta_1 )/({4 {\rm \pi}k_0^2 k^2}), \end{gather}$$

(3.7f)$$\begin{gather}\hat\varPhi_{23}^{Mann}(\boldsymbol{k}) = {E (k_0 )} ({-}k_2 k_{30}+ (k_1^2+k_2^2 ) \zeta_2 )/({4 {\rm \pi}k_0^2 k^2}), \end{gather}$$

where $k_0$ is the magnitude of the initial wavenumber vector $\boldsymbol {k_0}=(k_1,k_2,k_{30})$ with $k_{30}=k_3+\beta k_1$, and where $\beta$ is the non-dimensional eddy lifetime. Mann proposed an eddy lifetime model that has the correct asymptotic behaviour when compared with observations in wall-bounded flows. The non-dimensional eddy lifetime $\beta$ was expressed in terms of the hypergeometric function $_2F_1$. However, Wilson (Reference Wilson1998) reformulated this expression in terms of the incomplete beta function $B$:

(3.8)\begin{equation} \beta=\frac{\sqrt{3} \varGamma}{k \ell}\left[B_{1 /(1+k^2 \ell^2)}\left(\frac{1}{3}, \frac{5}{2}\right)\right]^{{-}1/2}, \end{equation}

which is much easier to compute using available numerical libraries.

The quantities $\zeta _1$ and $\zeta _2$ are given as

(3.9a,b)\begin{equation} \zeta_1=C_1-\frac{k_2}{k_1} C_2, \quad \zeta_2=\frac{k_2}{k_1} C_1+C_2, \end{equation}

with

(3.10)\begin{equation} C_1=\frac{\beta k_1^2(k_0^2-2 k_{30}^2+\beta k_1 k_{30})}{k^2(k_1^2+k_2^2)} \end{equation}

and

(3.11)\begin{equation} C_2=\frac{k_2 k_0^2}{(k_1^2+k_2^2)^{3 / 2}} \arctan \left[\frac{\beta k_1(k_1^2+k_2^2)^{1 / 2}}{k_0^2-\beta k_{30} k_1}\right]. \end{equation}

Unlike the HGW model, the Mann model is anisotropic, but it remains homogeneous in all directions. To account for the vertical inhomogeneity introduced by the wall, Mann (Reference Mann1994) proposed a modification to his original model based on a similar approach to that described in § 3.2. However, because of the non-zero mean shear in the model, the same procedure as in § 3.2 leads to a highly complex mathematical model that requires numerical integrations over the $k_3$ direction, with a computational cost of $\mathcal {O}(N_3^3)$ per horizontal wavenumber pair $(k_1, k_2)$. This would be a bottleneck in our assimilation algorithm. To avoid this complexity, we directly obtain ${\hat {R}_{ij}^{Mann}}(k_1, k_2, r_3)$ by applying the inverse fast Fourier transform to the model in (3.7). In contrast to the analytical IFT conducted over an infinite domain in the HGW model, the numerical IFT here yields a periodic correlation function in the vertical direction. This periodicity needs to be eliminated to avoid any unwanted effects on the assimilation results. To achieve this, the IFT is applied on an extended domain in the vertical direction, which is at least twice as large as the reconstruction domain height. The additional part is then discarded. A similar approach was adopted in Mann (Reference Mann1998) in the context of synthetic turbulence generation in a finite domain. While it is acknowledged that employing the homogeneous version of the Mann model, as opposed to the vertically inhomogeneous variant, may yield less precise results in close proximity to the wall, it should be noted that this discrepancy diminishes as one moves away from the wall. In this case, the two versions of the model are anticipated to exhibit similar behaviour, as indicated by Mann (Reference Mann1994).

Before utilizing the Mann model in our assimilation algorithm, one remaining issue needs to be addressed. Specifically, the model in (3.7) is undefined when $k_1=0$ due to the singularity of $\zeta _1$ and $\zeta _2$ in (3.9a,b). To avoid this issue, the limit values $\lim _{k_1 \rightarrow 0} \zeta _1=-\beta$ and $\lim _{k_1 \rightarrow 0} \zeta _2=0$ are used (Gilling Reference Gilling2009).

5.1. Set-up of lidar and virtual measurements

In this study, the measurement series ${\boldsymbol{\mathsf{Y}}} \triangleq [\boldsymbol {y}_0, \ldots, \boldsymbol {y}_{N_s}]$ is obtained virtually on a fine reference simulation (see § 5.2) using the implementation of Bauweraerts & Meyers (Reference Bauweraerts and Meyers2020) for the Lockheed Martin WindTracer lidar (Krishnamurthy et al. Reference Krishnamurthy, Choukulkar, Calhoun, Fine, Oliver and Barr2013). A full description of the implementation and set-up can be found in the original papers. However, for the sake of completeness, a brief overview is provided in this section.

In the current work, we collect our measurements in plan-position-indicator (PPI) scanning mode with a zero elevation angle and constant azimuthal sweeping covering the horizontal sector shown in figure 1. The measurements are collected for a period of $T_m=0.1 H/u_*$ from a single lidar mounted at $\boldsymbol {x}_m=[0,0,0.1H]$, with a range gate of $\Delta r=105$ m and total number of gates of $N_r=100$. Therefore, at each sampling time $t_n=t_0+nT_s$, with $T_s=5 \times 10^{-4}H/u_*=1$ s, a vector of wind-speed readings $\boldsymbol {y}_n=[y_{n,1},\ldots,y_{n,N_r}]$ is measured using ${y}_{n,i}={h}_{n,i}(\boldsymbol {u}({x}, t))$ and stored at equally spaced locations along the lidar beam. The observation function ${h}_{n,i}$ is given as (Bauweraerts & Meyers Reference Bauweraerts and Meyers2020)

(5.1)\begin{equation} h_{n, i}(\boldsymbol{u}(\boldsymbol{x}, t)) \triangleq \frac{1}{T_s} \int_{t_{n-1}}^{t_n} \int_{\varOmega} \boldsymbol{u}(\boldsymbol{x}, t) \boldsymbol{\cdot} \boldsymbol{e}_l(t) \mathcal{G}\left(\boldsymbol{Q}(t) \left(\boldsymbol{x}-\boldsymbol{x}_i(t)\right)\right)\, \mathrm{d}\kern0.7pt \boldsymbol{x}\, \mathrm{d}t, \end{equation}

where $\mathcal {G}$ is the lidar filter kernel, $\boldsymbol {e}_l$ is the lidar beam direction, $\boldsymbol {Q}$ is a rotation matrix from the reference to the lidar coordinate system and $\boldsymbol {x}_i$ is the measurement location given as $\boldsymbol {x}_i(t)=\boldsymbol {x}_m+(r_0+\Delta r(i-1)) \boldsymbol {e}_l(t)$, with the lidar range gate $\Delta r=105$ m and $i=1, \ldots, N_r$. The lidar beam moves in a horizontal sweeping pattern with a constant azimuthal angular velocity of $|\partial \phi /\partial t|=2\Delta \phi /T_m=0.00371$ rad s$^{-1}$, where $\Delta \phi$ is the azimuthal range covered by the lidar. Figure 1 summarizes this picture.

Figure 1.

Schematic of the reference (a) and reconstruction (b) domains with the lidar mount (purple) and the scanning region (inside dashed lines). The figure is reprinted from Bauweraerts & Meyers (Reference Bauweraerts and Meyers2020).

5.2. Set-up of reference simulation

In the current study, we assume a pressure-driven boundary layer with a height of $H=1$ km, a friction velocity of $u_*=0.5$ m s$^{-1}$ and a roughness length of $z_0=0.1$ m. After a spin-up period of $50H/u_{*}$ on a reference domain $30H \times 5.4H \times H$ with a fine grid of size $2000 \times 360 \times 200$, a series of virtual measurements are collected in the PPI scanning mode over a period of $T_m=200$ s and also $400$ s. After obtaining the measurement series ${\boldsymbol{\mathsf{Y}}}$, the reconstruction problem can be started. The reconstruction is done on a shorter domain with a coarser mesh than the reference simulation. In this study, the reconstruction domain has dimensions of $18H \times 5.4H \times H$ with a grid size of $360 \times 108 \times 60$ (see figure 1). Before solving the optimization problem (2.4), the background correlation tensor $\hat {\boldsymbol {R}}$ is constructed from the analytical models in § 3 on the reconstruction grid and stored to be used as a prior for the MAP problem in (2.4). Computing and storing the analytical tensor takes only a few seconds using the same computational nodes allocated for the optimization problem (see § 5.4). Therefore, they are usually deleted after solving the optimization problem and recomputed again when needed.

5.3. Set-up of correlation tensor

5.3.1. The HGW model

As discussed above, the HGW correlation tensor in (3.6) is a function of the length scale $\ell$ and the isotropic variance $\sigma _{iso}^2$. Meteorological data from the ABL (Kaimal et al. Reference Kaimal, Wyngaard, Izumi and Coté1972) show that the variances of the different velocity components $\sigma _{u}^2$, $\sigma _{v}^2$ and $\sigma _{w}^2$ in the surface layer (i.e. $x_3 \le 0.1H$) are not equal to each other (see table 1). However, because of the isotropic assumption in the HGW model, the model cannot pick up a preferred direction, meaning that only a single choice of the variance can be imposed. In this study, we choose the isotropic variance to be equal to the streamwise variance of Kaimal et al. (Reference Kaimal, Wyngaard, Izumi and Coté1972), leading to $\sigma _{iso}^2=4.77 u_*^2$. On the other hand, the length scale $\ell$ is chosen to match the inertial subrange asymptote (Peltier et al. Reference Peltier, Wyngaard, Khanna and Brasseur1996), which leads to the following expression (Wilson Reference Wilson1997):

(5.2)\begin{equation} \ell = 0.8743 \frac{\sigma_{iso}^3}{\epsilon}, \end{equation}

where $\epsilon$ is the turbulent kinetic energy dissipation rate. Using the similarity theory in the surface layer of a neutrally stable ABL, the dissipation rate can be approximated as (Stull Reference Stull1988)

(5.3)\begin{equation} \epsilon=1.24 \frac{u_*^3}{\bar{\kappa} z}, \end{equation}

where $\bar {\kappa } = 0.41$ and $z$ is the vertical distance, which is chosen to be equal to the lidar height $z=0.1H$. The resulting values of the parameters are summarized in table 3. It is worth mentioning that we noticed that the reconstruction results are not very sensitive to small changes in the model parameters.

Table 1.

Values of velocity variances as given in Kaimal et al. (Reference Kaimal, Wyngaard, Izumi and Coté1972).

As briefly mentioned in § 3, the analytical IFT for the spectral tensor ${\boldsymbol {\varPhi }}^{iso} (\boldsymbol {k})$ used in the HGW model is only available for the Batchelor turbulence case (see Appendix D), which is the main advantage of this model. Therefore, we limit our analysis for the HGW to the Batchelor case.

5.3.2. The Mann model

In addition to the two parameters discussed above ($\sigma _{iso},\ell$), the Mann model has an additional parameter $\varGamma$ which reflects the amount of distortion applied to the initial isotropic tensor. If $\varGamma =0$, no shear is imposed and the isotropic tensor in (3.3) is retrieved. As $\varGamma$ is increased, $\sigma _{u}^2$ and $\sigma _{v}^2$ increase while $\sigma _{w}^2$ and $\sigma _{uw}$ decrease. This behaviour is illustrated in figure 4 of Mann (Reference Mann1994). To choose the model parameters, Mann used a least-squares fitting algorithm to match the one-dimensional spectra produced by the model to the ones obtained experimentally. Here, we use a simpler procedure, similar to the one reported by Wilson (Reference Wilson1998), that is also more in line with the typical data that would be available from atmospheric measurements. Firstly, using figure 2, $\varGamma$ is chosen to satisfy the ratio $\sigma _u^2/\sigma _w^2$ obtained from the reference LES at $x_3=0.1H$ (see table 2). Afterwards, $\varGamma$ and $\sigma ^2_u/u_*^2$ are used to obtain $\sigma _{iso}^2$ from figure 2. Finally, the length scale $\ell$ is obtained similarly to the previous section. Table 3 summarizes the parameters obtained following this procedure for both Batchelor and Saffman energy spectra.

Table 2.

Values of velocity variances as calculated from the reference LES at $x_3=0.1H$.

Table 3.

Parameters of the spectral tensors in § 3.

Figure 2.

The normalized velocity variance of the Mann model as $\varGamma$ changes using Batchelor and Saffman energy spectra. The blue colour represents $\sigma _u^2/\sigma _{iso}^2$, while the red and green represent $\sigma _v^2/\sigma _{iso}^2$ and $\sigma _w^2/\sigma _{iso}^2$, respectively. The vertical dotted lines correspond to the selected $\varGamma$ values following the procedure in § 5.3.2 for both Batchelor and Saffman cases.

As will be described in § 5.4, the absolute value of $\sigma _{iso}^2$ in table 3 becomes somehow redundant in our 4D-Var problem due to the Pareto front technique used to fix $\gamma ^2$. Nonetheless, we opted to propose a way to fix it as part of the parameter selection process in the current section.

5.3.3. Two-point correlation tensor

In this section, we compute the diagonal component of the normalized two-point correlation function in the streamwise direction, given as $C_{11}(\boldsymbol {x},\boldsymbol {\breve {x}})= R_{11}(\boldsymbol {x},\boldsymbol {\breve {x}})/ \sigma ^2_{u}(\boldsymbol {x},\boldsymbol {\breve {x}})$, for both the Mann and the HGW models presented in § 3, and compare them with those obtained by Bauweraerts & Meyers (Reference Bauweraerts and Meyers2020) in figure 3.

Figure 3.

Cross-section of the two-point correlation $C_{11}(\boldsymbol {x},\boldsymbol {\breve {x}})$ with reference point $x = [0, 0, H/4]$ at $x_3 = H/4$ (b), $x_2 =0$ (a) and $x_1=0$ (c). The contour lines are drawn for $C_{11}=[-0.1,-0.05,0.05,0.1,0.3]$: solid black, $C>0$; dashed red, $C<0$. The figure is reprinted from Bauweraerts & Meyers (Reference Bauweraerts and Meyers2020).

Figures 4 and 6 present $C_{11}$ for the HGW and Mann models, respectively, for a reference point located at $\boldsymbol {x}=[0,0,H/4]$, which is a bit higher than the lidar height in § 5.1. However, we opted to show the correlation function at this height to facilitate comparisons with figure 3. Moreover, the difference between the two heights is negligible (not shown). The figures provide a visual comparison of the correlation functions between the two models. For the HGW model with Batchelor spectrum (figure 4), the correlation consists of three small positive central lobes with a spanwise width close to what is observed in Bauweraerts & Meyers (Reference Bauweraerts and Meyers2020). However, the HGW model underpredicts the correlation length by an order of magnitude when compared with the value of $10H$ reported by Fang & Porté-Agel (Reference Fang and Porté-Agel2015) around the reference point.

Figure 4.

Cross-section of the HGW two-point correlation $C_{11}(\boldsymbol {x},\boldsymbol {\breve {x}})$ with reference point $x = [0, 0, H/4]$ at $x_3 = H/4$ (b), $x_2 =0$ (a) and $x_1=0$ (c).

On the other hand, the Mann model in figures 5 and 6 exhibits a more complex correlation function. In the horizontal plane, the function displays three closed lobes elongated in the $x_1$ direction, with a size significantly larger than in the HGW case. However, it can be clearly seen that utilizing the Saffman spectrum (figure 6) leads to the longest correlations among the proposed models. In the spanwise direction, the correlation function decreases until it reaches a saturation point at a small negative value. In contrast, Bauweraerts & Meyers (Reference Bauweraerts and Meyers2020) observed long alternating positive and negative lobes extending throughout the entire domain in the streamwise direction, with a correlation width much less than the one observed here. These long lobes were attributed to the large-scale coherent structures in a neutral ABL as observed in Alcayaga et al. (Reference Alcayaga, Larsen, Kelly and Mann2020).

Figure 5.

Cross-section of the Mann (with Batchelor spectrum) two-point correlation $C_{11}(\boldsymbol {x},\boldsymbol {\breve {x}})$ with reference point $x = [0, 0, H/4]$ at $x_3 = H/4$ (b), $x_2 =0$ (a) and $x_1=0$ (c).

Figure 6.

Cross-section of the Mann (with Saffman spectrum) two-point correlation $C_{11}(\boldsymbol {x},\boldsymbol {\breve {x}})$ with reference point $x = [0, 0, H/4]$ at $x_3 = H/4$ (b), $x_2 =0$ (a) and $x_1=0$ (c).

In the $x_2\unicode{x2013}x_3$ direction, we also observe central positive lobes. In the $x_1\unicode{x2013}x_3$ plane, the Mann model produces positive lobes that are inclined in the direction of the mean flow with an inclination angle of around $19^\circ$, which is slightly higher than what is observed in other studies (Marusic & Heuer Reference Marusic and Heuer2007; Sillero, Jiménez & Moser Reference Sillero, Jiménez and Moser2014), where a value of around $10^\circ$ was reported. We observed that the inclination is very sensitive to the choice of $\varGamma$. However, in the current work, we prioritize a practicable way of selecting parameters, as outlined in § 5.3.2, rather than further fitting on the inclination angle, which is not directly identifiable as an input to the prior in our 4D-Var formulation.

5.4. Set-up of the assimilation problem

After determining the correlation tensor in the previous section, one remaining component to be chosen before solving the optimization problem (2.4) is the variance $\gamma ^2$. Looking at the optimization problem, the objective function can be rewritten in a more explicit way as $\mathscr {J}(\boldsymbol {\varTheta }_0)=\boldsymbol{\varTheta} _0^T \boldsymbol {B'}^{-1} \boldsymbol {\varTheta }_0/2\sigma _{iso}^2+ \sum _{n=1}^{N_s} \| \boldsymbol {y}_n-\boldsymbol {h}_n(\mathcal {M}_t(\mathcal {R}(\boldsymbol {\varTheta }_0)) \|^2 / 2 \gamma ^2$, where $\boldsymbol {B'}$ is the normalized correlation tensor. As can be noticed, the relative weight of each term with respect to the other is determined by the ratio $\gamma ^2/\sigma _{iso}^2$, instead of the individual values of $\gamma ^2$ and $\sigma _{iso}^2$. However, since $\sigma _{iso}^2$ was already fixed in § 5.3 as part of the parameter selection procedure, $\gamma ^2$ now serves as the only weighting parameter between the two terms in the cost function, and the absolute value of $\sigma ^2_{iso}$ becomes somehow redundant. For a fixed $\sigma ^2_{iso}$, as $\gamma ^2$ is increased, less trust is given to the likelihood term, leading to a smoother solution. Decreasing $\gamma ^2$ gives more significance to the likelihood term, and the obtained solution becomes more irregular. In order to choose a suitable value for $\gamma ^2$, a trade-off between solution complexity and smoothness needs to be achieved. Similarly to Bauweraerts & Meyers (Reference Bauweraerts and Meyers2020), this is done through Pareto front analysis, which requires solving the optimization problem multiple times per model for different values of $\gamma ^2$. Figure 7 shows that $\gamma ^2=10^{-3}$ achieves a good trade-off for both models, so it will be used in the rest of this study. According to our experience, the reconstruction accuracy appears to be much more sensitive to big changes in $\gamma ^2$ than the convergence behaviour.

Figure 7.

Pareto front for the optimization problem (2.4) using the HGW (a) and Mann regularization with Batchelor spectrum (b). The figure is obtained by solving the optimization problem many times using different values of $\gamma ^2$. The points are labelled with the corresponding $\gamma ^2$ value.

After choosing $\gamma ^2$, the optimization problem (2.4) is solved for an optimization time horizon equivalent to the lidar scanning period $T_m$, with an initial condition $\hat {{\boldsymbol {\varTheta }}}_0=0$. The spatial mean flow profile of the initial field as well as the surface roughness and $u_*$ are assumed to be known from the reference simulation. Starting from the initial condition, the optimization algorithm converges to a local minimum where $\boldsymbol {\nabla } \mathscr {J}=0$. In this study, the optimization algorithms stop upon achieving a relative value of the gradient norm $\| \boldsymbol {\nabla } \mathscr {J} \| / \| \boldsymbol {\nabla } \mathscr {J}_0 \|=6 \times 10^{-3}$, where $\boldsymbol {\nabla } \mathscr {J}_0$ is the initial gradient obtained using an initial guess for the initial condition $\hat {{\boldsymbol {\varTheta }}}_0=0$. Figures 8 and 9 show the convergence history of the HGW and Mann models with the Batchelor spectrum. Each iteration corresponds to an L-BFGS-B iteration, which involves a forward and adjoint simulation as well as a line search procedure (if the initial step $\alpha _k$ does not satisfy the Wolfe conditions). However, we only needed to perform the additional line search procedure a few times, so the total number of iterations mainly represents the number of forward simulations (or half of the forward and adjoint simulations). Figure 9 shows that we achieve a super-linear convergence rate when either of the analytical models is used.

Figure 8.

Cost function convergence (100 iterations) using the HGW (a) and Mann regularization with Batchelor spectrum (b).

Figure 9.

Convergence of the relative gradient using the HGW regularization (blue) and Mann regularization with Batchelor spectrum (orange).

As a point of comparison, we also investigated ridge regression. We found that it leads to slow convergence (a slope of $-1.1$), and only provides reasonable results inside the scanning area, as further discussed in Appendix E. We will use results from the ridge regression as a further point of comparison in the reconstruction error analysis discussed in § 6.2.

The optimization in this study was performed on the wICE supercomputer of the Flemish supercomputer centre using four computational nodes, each of which has two Intel Xeon Platinum 8360Y 2.4 GHz CPUs with 256 GB RAM. The total wall time of one assimilation problem with $T_m=0.1H/u_*$ (provided 100 iterations) is around 8 hours. It must be highlighted that this cost includes the preprocessing cost (involving the computation of the background tensor), which is negligible in our assimilation algorithm compared with the previous study of Bauweraerts & Meyers (Reference Bauweraerts and Meyers2020).