Low level winds weaker than U10 regionally in L58

[megandevlan]

At the CAM7 meeting, we discussed some new results comparing the magnitude of the lowest model level winds to the diagnosed U10 from the coupler. In a series of F2010climo cases, we found that in L32 and L48, regardless of dynamical core or physics package used, U10 is weaker than the lowest model level winds globally. But when moving to L58, that relationship can reverse:

The situation where U10 > windBot shouldn't be possible if the formulation of U10 is correct, but we see this consistently over the Southern Ocean in JJA. In DJF, the Southern Ocean still shows signs of the behavior but now the northern hemisphere is in play too; so ultimately, the winter hemisphere seems to be where the problem is primarily occurring.

Accompanying this is an increase in U10 and LHFLX biases, and an overall increase in the surface stress.

Per Thomas' suggestion, we're trying out a run with the COARE algorithm in shr_flux_mod rather than the default, so I should be able to add that update soon.

[megandevlan]

At least at first glance, using the COARE algorithm @tto061 suggested doesn't completely solve the issue of U10 > windBot in the Southern Ocean, though it does move things in the right direction.

Comparing the two rightmost columns in the fig below ('camdev_move' vs. 'camdev_move_COARE') illustrates the impact of switching from the default formulation of the ocean-atmosphere exchange code to the COARE version. There's still a surprisingly persistent negative sign over the Southern Ocean, though it's definitely weaker. The expected positive relationship between the two variables has also strengthened over the tropics (i.e., there's more of a positive difference between the 10m and ~22m winds when using COARE).

In terms of the U10 bias though, the COARE algorithm does perform better than the default option in most regions. One exception to that though is the NW Pacific, where biases have increased slightly.

Again the concern is that it shouldn't be possible for the lowest model level windspeed to be weaker than the 10 m winds, and the fact that this shows up across two algorithms is pretty surprising to me.

[JulioTBacmeister]

Does "camdevMove" refer to moving the output routine to something like its earlier position - before coupling?

[megandevlan]

Correct, "camdevMove" is basically identical to camdev, but I've moved the history write to just before coupling so that we can see the actual UBOT and VBOT that the coupler is using to compute U10 with.

[tto061]

Hi Meg could it be a time-averaging effect? u10 is computed as the time-mean of sqrt(duu10) where duu10=u10**2, so it is the true time mean of the wind speed. When you compute sqrt(ubot**2+vbot**2) you might be taking time-means for ubot and vbot. Is that the case? If so, the result will always be less than the time mean of the squares, and in areas of large wind variability the effect will be particularly pronounced and trump the (intrinsically) lower value of u10. Thomas Thomas Toniazzo Meteorologiska institutionen (MISU) Arrhenius Laboratory, Svante Arrhenius Väg 16C Stockholms universitet SE-106 91 Stockholm Sverige tel. +46 (0)8 16 2416

…

On 2022-05-03 16:13, Meg Fowler wrote: Correct, "camdevMove" is basically identical to camdev, but I've moved the history write to just before coupling so that we can see the actual UBOT and VBOT that the coupler is using to compute U10 with. — Reply to this email directly, view it on GitHub <#110 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADZGLJD7AKKABUYMWW6FQYLVIEX77ANCNFSM5U37P5ZA>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[JulioTBacmeister]

Hi Thomas, I have a Southern Ocean SCAM case forced by ERA-I derived tendencies which exhibits occurrences of u10 > sqrt( ubot^2+vbot^2 ) within the shr_flux_atmocn code itself (with the original iteration) , as well as in the instantaneous time-step output. So I think there really is something strange going on in the iteration.

[tto061]

Hi Julio! thanks, that makes it clear enough then. But let me try a last shot in this direction: are we sure about the sign conventions of us and vs, which are subtracted from ubot,vbot before computing u10? in case they're not zero, of course. Over the SO they (us,vs) should be noticeably non-zero and generally aligned with the (ubot,vbot) -- which should make u10 even smaller compared to |(ubot,vbot)|, unless the signs are wrong... I'm asking because funnily enough this did go wrong in some models in the past (i.e. the ocean current ended up being added to the wind rather than subtracted). But it's a long shot. I'll have another look at the iteration. Thomas Thomas Toniazzo Meteorologiska institutionen (MISU) Arrhenius Laboratory, Svante Arrhenius Väg 16C Stockholms universitet SE-106 91 Stockholm Sverige tel. +46 (0)8 16 2416

…

On 2022-05-03 17:54, JulioTBacmeister wrote: Hi Thomas, I have a Southern Ocean SCAM case forced by ERA-I derived tendencies which exhibits occurrences of u10 > sqrt( ubot*2+vbot*2 ) within the shr_flux_atmocn code itself (with the original iteration) , as well as in the instantaneous time-step output. So I think there really is something strange going on in the iteration. — Reply to this email directly, view it on GitHub <#110 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADZGLJFPCO7UWSE5GYDXRKDVIFD4NANCNFSM5U37P5ZA>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[swrneale]

We could test by commenting out the us,vs in a run. But thinking about it, these are AMIP runs so I am assuming there is no 'data model' for surface ocean drift, which I believe is correct. I think we need to dive into the TBOT, QBOT that is being passed into the shr_flux_mod routine, and examine the stable/unstable classification and whether the vertical stability functions are off the rails. I had set up an offline version of this code and I think Meg is going to just print out some instantaneous data from a short AMIP run to see how it behaves.

[JulioTBacmeister]

I agree that US and VS are likely zero for any FCASE runs. I'm pretty sure I've printed these out before and seen that this is the case. But I'll double check in my SCAM run.

[megandevlan]

To Thomas' point re: time averaging, it's possible there's some role of that, but as Julio mentioned it's still showing up even if we control for that. In the multi-year global cases shown above, I've been relying on daily mean U and V compared against daily mean U10. Right now though, I'm looking at some instantaneous output as Rich suggested, and I still see this behavior at various time steps (example below). As a side note, this plot should be compared to the DJF plot above not JJA (initialized runs on Jan 1st), so the signal we should see is more in the northern hemisphere vs the Southern Ocean.

I'm looking at the last 5 days from a week-long run branched from the global versions of these two resolutions for now, and I think the next step is getting this output run through Rich's offline script to diagnose the terms going into the stability function more directly.

[tto061]

OK, here I see a spatial match between areas of very low wind and negtive ubot-u10. This can be due to the 0.5 m/s "gustiness" added (quadratically) to the windspeed for the calculation of ustar and of duu10. So how would these plots look if you take sqrt(ubot**2+vbot**2+.25)? Thomas Toniazzo Meteorologiska institutionen (MISU) Arrhenius Laboratory, Svante Arrhenius Väg 16C Stockholms universitet SE-106 91 Stockholm Sverige tel. +46 (0)8 16 2416

…

On 2022-05-03 20:01, Meg Fowler wrote: To Thomas' point re: time averaging, it's possible there's /some/ role of that, but as Julio mentioned it's still showing up even if we control for that. In the multi-year global cases shown above, I've been relying on daily mean U and V compared against daily mean U10. Right now though, I'm looking at some instantaneous output as Rich suggested, and I still see this behavior at various time steps (example below). As a side note, this plot should be compared to the DJF plot above not JJA (initialized runs on Jan 1st), so the signal we should see is more in the northern hemisphere vs the Southern Ocean. Screen Shot 2022-05-03 at 11 49 40 AM <https://user-images.githubusercontent.com/42625059/166510303-0e4b4a5b-62a7-4668-9150-2e8e1c3898a1.png> I'm looking at the last 5 days from a week-long run branched from the global versions of these two resolutions for now, and I think the next step is getting this output run through Rich's offline script to diagnose the terms going into the stability function more directly. — Reply to this email directly, view it on GitHub <#110 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADZGLJFDYPQW5R5NJRRWYKDVIFSWVANCNFSM5U37P5ZA>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[megandevlan]

Ah, sorry for not clarifying this part earlier; in the two maps of instantaneous output, these are still using the default shr_flux_mod calculation, not the COARE algorithm. So in this case there isn't that gustiness term being added quadratically. I can re-run this with the COARE algorithm as well to confirm though.

@JulioTBacmeister, are your Southern Ocean SCAM runs that show this using COARE or the default formulation of surface fluxes (or is this signal coming through in both?)

[JulioTBacmeister]

My SCAM simulations are using the default surface flux calculation not COARE.

[megandevlan]

(Apologies for the lengthy update)
We have a mostly working copy of the offline shr_flux_mod calculations that Rich alluded to yesterday. I'll add a caveat though that the current version does not have the cold air outbreak option turned on (and it is on by default in the global runs), so these aren't quite a direct comparison. But the difference between U10 as computed by the coupler vs by this offline version aren't very different in most of the globe:

One point to be aware of though is that in the regions where we were seeing some of the biggest U10 > windBot signals is where the cold air outbreak addition to vmag is likely active (based on some plots of the criterion to activate it during these short simulations). I'll work on getting this functionality added next week so we can compare more clearly.

All this to say - when we're looking at U10 > windBot in these particular cases, we shouldn't focus on the largest signal in the Labrador Sea. But taking a mean over ~10 days in January, there are subtle signs of the signal over the northern Indian Ocean as well:

Looking through some of the terms that go into calculating U10 (again, just in the default version of the shr_flux_mod formulation), there's at least two interesting terms that pop out: "hol" and "rd_ov_rdn", with differences between the L58 and L32 resolutions in this Indian Ocean region. Hol is the height at ZBOT over L, with L being the Monin-Obukhov length scale. It would obviously have the potential change as the height of ZBOT's lowered, but it's also a pretty key quantity as it factors into the stability function that scales atmospheric winds down to 10m. The other term is the ratio of rd to rdn (momentum exchange coefficients). This ratio directly multiplies the wind magnitude vmag to give you u10n: u10n = vmag * rd/rdn. So if you increase this ratio to be >1 more often, you'll wind up with 10m wind speeds that are larger than the winds at the lowest model level.

So ultimately the patterns in hol and rd/rdn differences tends to line up closely with where UBOT>windBot, if we exclude where the cold air outbreak functionality is likely at play. It's nothing terribly concrete just yet, but at least gives us an idea of where to dig further. It's also a bit suspicious that the largest signals are where the cold air outbreak is getting turned on; I don't know if that's going to be the case in the longer runs as they have a much larger region of U10>windBot anyway, but it's worth looking into more (maybe simulate a few years with the option turned off).

[swrneale]

Hi Meg, this all makes sense. I'm not sure how to calculate this 10/L term though as it seems a T and Q is required for the ZBOT/L term which we have, but what do we use for 10/L? One further thing that occurred to me is that u10n in the model is a neutral layer windspeed so we may not necessarily see the UBOT and u10 values converge as the lowest model level gets near 10m. That being said that, what happens if you redo your figures above and set zbot=15 or even 11m, do the functions behave?

[adamrher]

Maybe it's time to reach out to Bill Large?

[JulioTBacmeister]

Does u10>ubot occur within shr_flux_atmocn when the COARE subroutine is used and no cold_air_outbreak enhancement occurs?

[megandevlan]

@swrneale Good point about the difficulty in computing the 10/L term. I'm assuming it's too simplistic to say that we have the value of ZBOT/L and can back out L from that? (Not sure how height dependent L is on the T and Q used to compute it, but if it's truly some estimate of the height at which buoyant production of turbulence equals shear production of turbulence, it should be robust).
Also correct that u10n is the neutral layer windspeed, though given how well U10 from the model agrees with the offline u10n computed for most of the globe I'd be surprised if this is a big problem. It looks like in med_phases_aofluxes_mod, u10 is actually directly from duu10n, which is just u10n^2 (though I could be missing a step somewhere):
aoflux_out%u10(n) = sqrt(aoflux_out%duu10n(n))
We could test just lowering zbot manually in these offline calculations though any BOT variables would then be off (TBOT, UBOT, etc. are still at 22m). Could be worth exploring though.

@JulioTBacmeister I'll have to code up the COARE subroutine in an offline form to find out, but I'll get started on that. I can say that at least within the default, u10>ubot is happening at various time steps within shr_flux_mod and in the mean sense over the Indian Ocean. Checking what happens with COARE is probably a good step to complete, and then we may want to think about @adamrher's suggestion of reaching out to Bill.

[JulioTBacmeister]

In addition to reaching out to Bill Large I suggest we all pore over the CAM4 documentation https://www.cesm.ucar.edu/models/ccsm4.0/cam/docs/description/cam4_desc.pdf
pages 157-162 (especially 161-162). This description is making a little more sense comparing with the code than Large&Pond 1981

[megandevlan]

It’s taking a bit longer than expected to code up an offline version of the COARE algorithm, but things are almost working in the tropics. In the extratropics, it looks like there’s some issues with the offline-computed u10 not quite agreeing with the U10 that’s output in the history files, with some odd banding behavior in the pattern of disagreement (difference averaged over the last 19 days of a 21 day run):

I’d also note though that this run is branched from year four of a simulation that applied the COARE algorithm, but the simulation used camdev physics and SE; the ones above (the L32 and L58 cases used to force an offline shr_flux_mod) used cam6 physics and FV, so not quite an even comparison. I wouldn't have thought that would contribute to this disagreement, but I'm not sure just yet.

In the instantaneous output written to the history files, I still find occurrences of windBot < U10 at various time steps, but the nature of that occurrence seems different than what we saw in the L58_fv_cam6 case:

The negative points here with COARE look like they occur along frontal boundaries, while the L58_fv_cam6 case showed a more widespread issue in certain regions (i.e. the North Indian Ocean) even at individual timesteps.

Unfortunately, the places where we see this negative signal in COARE is where the offline computed sqrt(duu10n) doesn’t match up with the U10 that’s output from the run, so I’d like to debug the offline code a bit more before I can say for sure that using COARE could be one path to eliminating the issue. One point of optimism for that possibility though is that if I take the instantaneous output from the L58_fv_cam6 simulation and run it through the COARE algorithm, we do seem to eliminate instances of u10n > vmag (or in COARE variables, sqrt(duu10n) > du). That assumes I’ve coded up the algorithm properly (still a little in question), but so far the signs are moderately encouraging.

[megandevlan]

Thanks Julio, this is a really nice way to visualize the issue. Over oceans in the L58 grid, zbot does tend to be near 20 m from what I can tell, so that's probably our biggest issue. There's a line in the 1982 Large and Pond explanation of that V(z) formulation (right after their Eq. 5) that I'm wondering about in light of this. After outlining that equation, they state "except at small Z when the logarithms of (z+z0)/z0 ...must be used." There's no definition of what a 'small' z value would be (and I would've assumed it to be on the order of a few meters at most), but perhaps we're beginning to bump into something similar? It would be nice to find mention of this somewhere else as well though, hopefully with a fuller explanation so I'll take a look around.

I'm not entirely sure I understand the LPY formulation either, but I think the roughness length is there (albeit hidden) in the calculation of rdn. Looking at the way z0 is defined in the CAM4 doc (thanks for suggesting we start digging through that, by the way):

If you were to plug that top equation into log(z/z0) you'd find something similar to what's in the current code (assuming I'm remembering my log rules correctly). The factor of 10 out front of the exp() term is essentially how we wind up with log(z/zref), and the function inside the exp() term is -k/cdn. I'm not entirely clear how we go from that to the U10 formulation, but at least there's decent similarity.

[swrneale]

@megandevlan so it seems the hol calculation needs to be broken open a bit. It involves tstar/tbot/thbot/qstar/ustar so would you plot these for L32/L58 offline to see how they may differ feeding into hol?

[megandevlan]

@swrneale Absolutely! Zooming into the main region with the signal, the fig below shows the difference in lowest level wind speed vs. the 10m wind speed computed in the offline shr_flux_mod (so vmag vs. u10n) on top for both the L32 and L58 cases. Where we see that signal become much more negative is where the hol value has become substantially more negative as well (middle three panels), in line with Julio’s plot and what we saw before:

If we want to break up the terms in hol, we need to consider as you say the values of THBOT (added code to output this directly from the model), QBOT, tstar, qstar, and ustar:
hol = loc_karman*loc_g*zbot(n)*(tstar/thbot(n)+qstar/(1.0_R8/loc_zvir+qbot(n)))/ustar**2

I’ve broken down those terms below with various degrees of granularity below. Ultimately, we have a substantial decrease in zbot (on the order of ~45 m in the North Indian Ocean region here), but that's pretty consistent spatially and thus can’t explain the regional variations in how much more negative hol becomes. It seems like that’s a combination of things: ustar^2 becomes smaller in this region so that the magnitude of hol can increase (become more negative), and there are some mixed changes in the numerator term as well that don't seem to be offsetting the decreased u*^2. I'd like to dig into these a bit more though, just initial thoughts so far.

Another important point that’s come up here before (Rich, I think you’ve pointed this out in particular) and that Julio’s reiterated recently: as far as we can tell, the value of U10 that’s being written out in the LPY algorithm is based on u10n, which is the the neutral equivalent wind. In unstable conditions, that value can be larger than the actual wind speed that would be expected at that height. It doesn’t look like the stability-corrected 10m wind is computed/used, so this could be complicating the comparison of windBot and U10. I'm not sure if the COARE algorithm also computes a neutral equivalent windspeed or if this is the actual stability adjusted value, but if the latter that would help explain why we don't really see 10m winds > windBot in that case.

[swrneale]

@megandevlan I think you are right about the neutral u10m being output, but that is the same for L32 and L58. And wouldn't the unstable u10m be even greater than this? Are the above figures base on the mean of several months or something? I wonder if you could plot the PDF at the center of that hol minimum to see if it is being affected by just a few hol values of -10 (the limit0 or something. In the end it is hard to see that the differences in hol between L32 and L58 can be ascribed to the component terms in hol. Keep going! :) !

[megandevlan]

It would be an issue with the L32 case as well, certainly, though I don’t think it would show up as strongly. Most of what we’re looking at is the difference between u10n and vmag, but the winds near 60m are of course stronger than the ones at ~20m in the L58 case. So I think the higher resolution puts us in a situation where u10n can exceed vmag much more frequently. That's my first thought, at least.

The above figures are based on about 12 days of a branched simulation in January, so it’s not a ton of times being averaged together. I agree though, these aren’t helping me make sense of the ‘hol’ terms either. I’ve been looking at a single gridcell where the vmag < u10n signal is pretty strong, and so also where the ‘hol’ decrease is pretty sizable:

At this particular point, even the L32 case can show u10n>vmag, but it’s a constant feature in L58. The limiter is absolutely coming into play more often with the increased vertical resolution, and that could be contributing to part of the signal. In general though, the L58 hol values also tend to be more negative than the L32 case, even if not hitting the limiter.

It's a bit simpler to see the contribution of the terms (though still not entirely straightforward) this way. The friction velocity (u*) has decreased substantially with the move to L58, pretty much across all time steps. In the code, this is updated in the iteration loop as ustar = rd * vmag; vmag has been reduced, but given that rd is shifted to the measurement height/adjusted for stability, I would’ve thought that would balance accordingly. Regardless though, the ustar^2 decline is more of a steady signal than the change in hol, so that’s not all that’s happening.

The humidity term (qstar/(1/zvir + qbot)) has increased its magnitude for the bulk of the period as well, while the temperature term (tstar/thbot) has become more negative only part of the time but notably this is when hol decreases most strongly (given the difference in magnitude between the temperature and humidity terms, this probably makes sense). So I think there’s a few underlying changes that are pretty constant and could driver larger negative hol values, but the signal overall seems pretty closely linked to the temperature term becoming more negative during distinct periods, forcing hol to bump into the -10 limiter more often.

We can play around with the limiter a bit in the offline function to test how much of an impact that has on this issue (probably worthwhile), though given that the times we hit -10 don’t necessarily seem to be a unique driver of u10n > vmag at this gridcell, I’m a little doubtful for now that this is the biggest issue we’re facing. I think as Julio’s plot shows, it’s stemming in large part from log(zbot/10) being so small now, too. But overall, next we can test the limiter’s role and we can look more into the impact of using u10n vs the stability adjusted u10.