The revised version of the paper is much improved. There is now a good motivation, the method description is stringent and the examples shown for the application are well chosen. However, in the method comparison S3D is not correctly applied. That is detailed in the major comment below and results for a correct application shown in the attached Figure. In addition, there are some inaccuracies in the description of the other methods. You should also point out the basic assumption of attributing all variance to a single wave together with the method overview. Some additional specific comments are given below. If the method intercomparison will be corrected and the other comments taken into account I recommend publication of the manuscript in ACP.
Major comment:
You are using an 1D data set. At this you cannot apply S3D or 3DST. There are two ways attacking the problem:
a) You can expand the data set to a 3D data set.
b) You can use an 1D sinusoidal fit in a sliding window
These two possibilities are discussed here for S3D, code for reproduction is attached and can be run with the general IDL package of the Juelich remote sensing group.
a) In order to produce a data set comparable to normal S3D application, expansion of the test function to 3D should include a wave structure in the two other spatial dimensions as well. We can choose two ways of doing this:
a.1) Use the same wavelength for both wave packets
a.2) Use different wavelengths for bot wave packets
Employing a cube-volume of the size of 0.5 in the direction of the variation, a.1 (red) is very close to the solution of the Hilbert transform, a.2 (blue) can separate between the wavelengths and has a step-wise transition in wavenumber without the spurious wavenumbers between the two solutions. The results are shown in the attached figure, panels a (amplitude) and b (wavenumber). However, for real-world problems we probably would apply a size of 1.0 or 1.5 in order to get more confidence on not ideally-shaped wave patterns. In the case of box-size 1.5 (panels c and d) we find a notable smearing of the envelope (15% underestimation of the peak-amplitude), still closer to the reference than the Stockwell transform. All these fits draw information from the two additional directions which we have chosen to be homogeneous in amplitude and wavelength.
b) A mere 1D sinusoidal fit of window length 1.5 produces somewhat smeared amplitudes and some spurious oscillations (panels e and f). That can also be emulated by S3D and prescribing no wavelength information in the two additional directions (A.3).
For the paper you could follow one of the following two approaches:
1. Use the 1D solution. Do not call it S3D. Mention that the additional dimensions will add to the quality of the results
2. Discuss a.1/a.2 for a realistic window length of 1.5. Note in the text that the window length has an influence.
Anyway, none of this approaches performs as lousy as the so-called S3D solution in the method intercomparison.
Similarly, the Stockwell transform could perform better than in this test. There is a tuning factor c to reduce the width of the Gaussian envelope and also for the ST the third dimension adds information. One could thus also tune ST for a closer match in the 3D case. Also this should shine up in the discussion.
Summarizing: the idea of a simple method intercomparison adds value to the paper, but it must be done right!
Minor comments:
P4L7 The importance of oblique propagation is not restricted to the mesosphere. Ehard et al show this for the mid stratosphere but it may be quite substantial already in the UTLS (e.g. Krisch et al). Oblique propagation hence just redistributes selective transmission and makes the analysis more complicated. It is definitely a limiting factor for the interpretation of vertical profiles and some discussion should be added.
P3L4 In principle you should be able to calculate time-space spectra of u and w and calculate cospectra for the respective wave modes (e.g. Alexander et al. 2004). This is not limited by the size of the analysis volume and can be used to apply zonal mean quantities. Given that one reviewer asked for the motivation to use horizontal divergence, that is still not sufficiently motivated. One motivation could be the correction terms which need to be applied for the influence of Coriolis force on pseudomomentum flux.
Figure 5c and interpretation: The vertical wavelength is shortest in the wind maximum and then increases(!) towards where you assume the critical level is. There is only a very small peak at that altitude (much longer wavelengths than at the assumed excitation altitudes), comparable in size to other local variations. There are a number of explanations beside critical level which could explain your finding of a maximum in the jet, for instance oblique propagation, partial reflection ... So, if you want to retain this discussion, you should offer more proof and infer the phase speed for the considered waves considering the phase shift between two consecutive time steps of ECMWF.
Specific comments:
P1L25 S3D without hyphen
P2L1 Ern and Preusse, 2012 and Ern et al. 2014 use different methods -> omit
P2L2 to P2L18
What you say is all correct, but I think you are missing the chance to set this in a broader context. The fundamental problem of all wave analyses methods is given by the fact that you have to compromise between spatial and spectral resolution (cf. the uncertainty relation). A second point is to which extent the sampling and total volume will influence the result. Going along such lines you can order along the spatial resolution:
Fourier Transform: Assumes homogeneity of amplitudes and phases through hole considered volume. Wavelengths determined by volume size. Covers whole spectrum continously. No information loss.
3DST: Assumes homogeneity inside a volume corresponding to one wavelength (or fraction, cf. scaling factor). Wavelengths limited by size of total analysis volume. Selecting largest events implies loss of variance (information).
S3D: Assumes homogeneity inside an predefined analysis volume. Restriction to wavelength by cube size (sensitivity study: reliable for wavelength < 2.5*cube size). Outcome depends on pre-selected cube size. Selecting leading (not largest) components implies loss of variance (information).
Hilbert Transform: Does not assume homogeneity. No limitation on wavelengths except Nyquist. Gives information of amplitude and phase, wavelength information inferred from local phase gradient -> all variance is attributed to one wave mode.
This prize you have to pay is important and should be mentioned for a fair assessment!
Anyway, if you present these different compromises between wave resolution and spatial resolution you can claim that Hilbert transform provides the one we are missing. UWADI thus provides the user with a complementary tool for investigating wave events.
All localized methods can be shifted by steps smaller than the analysis method, so one is, in principle, able to calculate wave parameters for each original point for 3DST and S3D. A large asset of the Hilbert transform is that it is computationally cheap.
More specific:
P2L2 "two" is most often chosen, but not a general limitation. However, as you say in your next sentence, always a small number is taken. I would just omit the "two".
P2L8 3DST, please use notations as in the reference papers
P2L11 The only assumption made is that of upward propagating waves. This ambiguity between e.g. eastward-upward and westward-downward cannot be decided from temperature observations alone and is not a restriction of the method.
P2L12 S3D uses the largest described variance, i.e. minimum chi-square to pick the wave, 3DST the largest amplitudes. In both cases it can be more than one wave. Anyway, that should not be relevant for what you want to say, perhaps:
Both S3D and 3DST use a small number of the most-prominent waves. This leaves some variance unattributed and hence means a loss of information.
What perhaps is more important in this context:
Both methods implicitly assume homogeneity of amplitude and wave vector, S3D inside one analysis volume, 3DST inside a volume corresponding to one wavelength.
This is a loss of spatial information.
P2L27 to make clear it is not flow over convection: by flow over orography, by convection ...
P4L18 use e.g. for the reference
P4L22 The phase of the wave ...
P6L4 Alienation of outliers is taken care of by two different quality. Please reformulate
P8L11 You can focus on these scales. However, a vertical limit of 15km will lead to shift part of the spectrum in and out the so-chosen visibility filter (Alexander, JGR, 1998 and later literature on this topic).
P8L27 do not interact with the mean flow. *W*ave action
It should not change for refraction either, please be more specific about the mode of interaction (e.g. breaking, dissipation, non-linear ...)
P9L8 In order to make the paper easier to read you could indicate by the notation whether a quantity is intrinsic or ground based; e.g. use \hat{\omega} and c_{g}.
P9L9 The wind also needs to be constant in time. Perhaps better turn around: If the wind ... Quote a ray-tracing paper, e.g. Lighthill or Marks and Eckermann.
P9L14 swap the two sentences: The displaced ... The jet streak ...
P9L16 I am wondering: is eastern Siberia really a comparable situation to Europe? The GW fronts seem more aligned than across the winds and there is less deceleration but strong curvature.
P9L23 However, it is counter-intuitive. Both a stronger damping and a coarser vertical grid would let expect the shorter waves to be more strongly damped. Did you try what happens if you increase the upper limit of your analysis range?
Fig 2: Omit ** on the exponents.
Fig 3: Please generate a second panel with zonal mean wind and N. At constant c_g, k vertical wavelength is inversely proportional to N, so you should see that as well.
Fig 4: It is evident, still: Give the units of the shown quantities, please.
P9L30 ? height range ?
P10L1 descented -> descended ?
It would be nice to indicate the profile locations in Fig 2a.
P10L7/L9 2* Overall
P10L10 ... remains constant. Puzzling, that's in the sponge.
P12L7 This could be one of the major assets of a local method - that one sees the scaling of the wave properties with the varying wind ...
P12L13 At these altitudes, GWs of vertical wavelength ...
P12L19 an SSW (i.e. an EsEsDoubleU)
P12L21 They point out theoretically that during the upward propagation of GWs these waves (comma, upward)
P12L22 Actually that was profile-based MEM/HA (Preusse et al, JGR, 2002) with a 10km running window.
For current day limb sounders we get only an estimate of the absolute horizontal wavelength and do not estimate wave action.
P12L25 the longest vertical wave length with 7 km
If that were really the longest existing wavelengths AIRS should not see any GWs at all. There should be also cases of substantially longer wavelengths.
P12L26 Scandinavian
P14L12 put -> exerted
P14L17 jet-generated GWs tend to have such phase speeds.
(convection would be much faster)
P14L19 relatively high
P14L27 Low wind speed and low wave action, maybe. In cases of larger waves these may still dominate the statistical effects leading to the higher wavenumbers.
P14L32 such -> thus ? or omit it
P15L2 nomination -> assumption or attribution of the variance to
Program code:
pro s3d_uwadi_test_3d,nzf=nzf
;
; Generates the UWADI-paper test function, adds a wave in x and y and
; runs an S3D
;
; Keywords:
; nzf number of points for the height of the S3D cube
loadct,39
if not keyword_set(nzf) then nzf= 15
divfa= 1.0 ; use 1.0 for box length of 1.5 and 3.0 for box length of 0.5
nz= divfa*fix(40*!pi)
z= dindgen(nz+1)/(nz*1.0)*4*!pi
x= dindgen(11) & y= dindgen(11)
setz_3d_field,x,y,z,xx,yy,zz,/double
; a.1 use the same wavevector for both wave packets
v= exp(-(zz-4.5)^2)*cos(4*zz+1.*xx+1.*yy)+exp(-(zz-7.5)^2)*cos(9*zz+1.*xx+1.*yy)
; a.2 use different wavevectors for both wave packets
; v= exp(-(zz-4.5)^2)*cos(4*zz+1.*xx+1.*yy)+exp(-(zz-7.5)^2)*cos(9*zz+1.*xx-1.*yy)
; emulate 1D, a.3
; v= exp(-(zz-4.5)^2)*cos(4*zz)+exp(-(zz-7.5)^2)*cos(9*zz)
a= exp(-(z-4.5)^2)+exp(-(z-7.5)^2)
plot,z,v[0,0,*], $
charsize=1.5,ytitle='value / amplitude', $
yrange=[-1,1]*1.0 ; 1.4
oplot,z,a,thick=2
hifiw= nzf/2
ampl= fltarr(nz+1)+!values.f_nan & chiq= fltarr(nz+1)+!values.f_nan & k= fltarr(3,nz+1)+!values.f_nan
; perform fits in a window/box running in z direction
for i= hifiw+1,nz-hifiw-1 do begin
xvx= reform(xx[*,*,i-hifiw:i+hifiw])
xv= fltarr(3,n_elements(xvx))
vv= fltarr(n_elements(xvx))
xv[0,*]= reform(xx[*,*,i-hifiw:i+hifiw])
xv[1,*]= reform(yy[*,*,i-hifiw:i+hifiw])
xv[2,*]= reform(zz[*,*,i-hifiw:i+hifiw])
vv[*]= reform(v[*,*,i-hifiw:i+hifiw])
print,i
; Use nested grids for finding best k
; sfk= sinfit_k_3d(reform(v[*,*,i-hifiw:i+hifiw]),xv,dx=[1,1,z[1]-z[0]],fiste=reform(v[*,*,i-hifiw:i+hifiw]))
; Use
sfk= sinfit_k_grad_3d(vv,xv,dx=[1,1,z[1]-z[0]],fiste=reform(v[*,*,i-hifiw:i+hifiw]))
; Fit with know k vector
; sfk= sinfit_raw_3d(vv,xv,[1.,1.,4.])
if (sfk.k[2] lt 0) then sfk.k*= -1.
ampl[i]= sfk.ampl & chiq[i]= sfk.chiq & k[*,i]= sfk.k
endfor
oplot,z,-a,thick=3
oplot,z,ampl,color=254,thick=3
oplot,z,-ampl,color=254,thick=3
stop
good= where(a gt 0.05)
plot,z,k[2,*],yrange=[0,20], $
charsize=1.5,charthick=3, $
ytitle='wavenumber'
oplot,z[good],k[2,good],thick=5
stop
return
end
pro s3d_uwadi_test,win_length=win_length,resolution=resolution
;
; Generates the UWADI-paper test function and runs a sinfit
;
; Keywords:
; win_length fit window length
; resolution dk/k when nested intervals are stopped
if not keyword_set(win_length) then win_length= 1.5 ; about 2*!pi/4.
nx= 1000
; nx= fix(40*!pi)
x= findgen(nx+1)/(nx*1.0)*4*!pi
v= exp(-(x-4.5)^2)*cos(4*x)+exp(-(x-7.5)^2)*cos(9*x)
a= exp(-(x-4.5)^2)+exp(-(x-7.5)^2)
loadct,39
plot,x,v, $
charsize=1.5,ytitle='value / amplitude', $
yrange=[-1,1]*1.0 ; 1.4
oplot,x,a,thick=2
oplot,x,-a,thick=2
stop
iwin_length= fix((win_length/(x[nx]-x[0])*nx)/2)*2+1 ; use odd number
hifiw= iwin_length/2
ampl= fltarr(nx+1)+!values.f_nan & lx= fltarr(nx+1)+!values.f_nan
for i= hifiw+1,nx-hifiw-1 do begin
sfk= sinfit_k(v[i-hifiw:i+hifiw],x[i-hifiw:i+hifiw],resolution=resolution)
ampl[i]= sfk.ampl & lx[i]= sfk.lx
endfor
oplot,x,ampl,color=254,thick=2
oplot,x,-ampl,color=254,thick=2
stop
good= where(a gt 0.05)
plot,x,2*!pi/lx,yrange=[0,20], $
charsize=1.5,charthick=3, $
ytitle='wavenumber'
oplot,x[good],2*!pi/lx[good],thick=5
stop
return
end
function sinfit_k,v,x,k_range=k_range,denser=denser,resolution=resolution
;
; Uses nested intervals to determine k in addition to amplitude and
; phase
;
; Keywords:
; k_range first guess for the interval, if not set use
; 2*!pi/[2.5*(range of x values),Nyquist limit]
; denser allow for denser sampling at first guess
; finer split each interval in finer new ones
if not keyword_set(finer) then finer= 10
if not keyword_set(resolution) then resolution= 0.0001
if not keyword_set(denser) then denser= 1.0
if not keyword_set(k_range) then begin
lenx= max(x)-min(x)
dx= lenx/(n_elements(x)-1)
k_range= 2*!pi/[2.5*lenx,2*dx]
endif
finer= finer*1.0 & denser= denser*1.0
; use about twice as dense as a Fourier-grid for the first try
dk= (k_range[1]-k_range[0])/(n_elements(x)*denser)
n_try= n_elements(x)*denser
k_try= dindgen(n_try+1)*dk+k_range[0]
k_best_idx= sinfit_k_chi_eval(v,x,k_try)
while (dk/k_try[k_best_idx] gt resolution) do begin
k_try_new= 0
if (k_best_idx eq 0) then $
k_try_new= dindgen(finer+1)*(k_try[1]-k_try[0])/finer+k_try[0]
if (k_best_idx eq n_try) then $
k_try_new= dindgen(finer+1)*(k_try[n_try]-k_try[n_try-1])/finer+k_try[n_try-1]
if (n_elements(k_try_new) eq 1) then $
k_try_new= dindgen(finer+1)*(k_try[k_best_idx+1]-k_try[k_best_idx-1])/finer+k_try[k_best_idx-1]
k_try= k_try_new & n_try= n_elements(k_try)-1
k_best_idx= sinfit_k_chi_eval(v,x,k_try)
dk= (k_try[1]-k_try[0])/finer
endwhile
return,sinfit(v,x,lx,/cal,k=k_try[k_best_idx])
end |
- Printer-friendly version
- Supplement
