Draxler Turbulence Model

[rudymcintyre]

Description

• Introduces TurbulenceModel base class to define difference turbulence behaviours.
• Defines the existing (default) turbulence implementation as AngularModel.
• Defines turbulence model based off Draxler calculations (DOI included in docstring)
• Takes sigma_hor/sigma_vert calculation out of compute_coupling_array during the forward model calculation and calculates the plume stability from the turbulence model in use.

Issues:

• I am not sure on the fact that wind_turbulence_horizontal now can represent values in either deg or m/s, but did it this way to avoid saving turbulence_deg and turbulence_ms in separate arrays This is arguably now a breaking change since the API of Meteorology has changed.
• meteorology.py still needs the function to calculate the turbulence in m/s. We calculate from high res wind and store as low res. Not sure how the use cases of this library would dictate the implementation of this.
• (22/04/2025): gaussian_plume.py still passes round one array for each of the wind turbulence parameters, giving us one array representing two units again. After the interpolation step, the turbulence is passed around as a arrays instead of inside the met object, so we'd be adding two more parameters to all of the 'inner' GP calculation functions. I will wait to see what the consensus is on this before proceeding.

Type of change

Please delete options that are not relevant.

• Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected). May cause conflicts with existing (external) code due to the change in name of the wind turbulence parameter arrays.
• This change requires a documentation update

Jupyter Notebooks

This model can be run with the basic Jupyter notebook included in the repo with the following changes:

Meteorology definition:

# Comment one or both of the flat value turbulence assignments. Adding Draxler vertical
# turbulence will require a 'w_component' array to be set though.

# met_object.wind_turbulence_horizontal_deg = 5.0 * np.ones_like(met_object.wind_direction)
met_object.wind_turbulence_vertical_deg = 5.0 * np.ones_like(met_object.wind_direction)

# Since the time axis is at a frequency of 120s, a large window is given here to include more data points.
met_object.calculate_wind_turbulence_horizontal_meter_per_sec(window="10min")

Gaussian Plume definition:

# Set the horizontal model (and/or vertical, depending on the content of your meteorology object)
# to be a DraxlerModel, using the appropriate default params.

# dispersion_model = GaussianPlume...
dispersion_model.turbulence_model_horizontal = DraxlerModel.default_horizontal()

How Has This Been Tested?

• Happy to write unit tests if you think this is a valuable contribution
• (22/04/2025). With the changes in gaussian_plume.py, there is definitely some more testable code now. I will add tests if you think this is the correct approach!

Checklist:

• My code follows the style guidelines of this project
• I have performed a self-review of my code
• I have made corresponding changes to the documentation
• My changes generate no new warnings
• I have added tests that prove my fix is effective or that my feature works
• New and existing unit tests pass locally with my changes
• Any dependent changes have been merged and published in downstream modules

[TannazH]

I suggest having the wind_turbulence_horizontal_deg, wind_turbulence_vertical_deg, wind_turbulence_horizontal_meter_per_sec (wind_turbulence_u_component) and wind_turbulence_vertical_meter_per_sec (wind_turbulence_v_component) as the properties of the meteorology object. I would also suggest implementing a method in the meteorology object to calculate the horizontal and vertical turbulences in m/s (similar to calculate_wind_turbulence_horizontal). I think this will make it clearer and we can use either option based on the need.

[rudymcintyre]

I have made this update in meteorology.py. There are now a lot of places where references should be updated, I will go through and update these but thought I'd push my changes so far anyway.

I have also implemented the horizontal and vertical wind turbulence parameter calculation functions. The perceived benefit of this turbulence model is a better estimate of turbulence using the HF wind (though we have seen similar results at 1Hz and 10Hz with this model). Since the current implementation of Meteorology only has one time axis, I propose that to use HF wind in this way, you could create an HF meteorology object and calculate the parameters there, then resample and assign to the LF met object.

I realise that this is a bit of an awkward workflow, but I think that perhaps the HF use case is quite niche, and so we shouldn't clutter the meteorology class with dual time axis artifacts. I would be interested to hear your thoughts.

[TannazH]

I was looking at the Draxler paper where you take the default parameters. I suggest referring to the table or the equation in which you obtained these values in the documentation.
I think it would be also beneficial to explain what each parameter means.
If one wants to override these default parameters, what would be the source to obtain them?
How sensitive are the estimated sigma_hor and sigma_ver to these input parameters in case one wants to use these default values instead of those dependent on the location and the ground threshold?

[rudymcintyre]

Looking back at this, a lot was really unclear. This is definitely not my specialist subject and Jason has not been around today to fill in any gaps but I've made an attempt at providing clarity based on what I can understand of the Draxler paper. I would appreciate your feedback.

In regards to your final points, I'm not sure I could recommend another source of parameters to use. There is this paper which we have also used for reference material, but seems to have less thoroughly defined parameters for the universal functions. (pp27)

Again, I am unsure about the sensitivity of the parameters to the sigma values, but I have attempted to provide some explanation to this effect, though it is not particularly rigorous.

[TannazH]

I suggest adding a minimal, reproducible code to the pull request as an example.

[rudymcintyre]

I have included the necessary changes to the example notebook in the PR description.

[bvandekerkhof]

Hi Rudy,

Given your comment regarding gaussian_plume.py about passing arrays which could be different units I feel it might be a bit confusing to upate all the inner plume functions with more potential parameters which specify if it's degrees or m/s, leaving it like it is now is also a bit strange I think as it feels like mix of notation. So, I think the implementation you had before where the parameter was named the same only the unit might be different was a better way after all. Would it be possible to revert back to that commit, hence not have meter_per_sec or deg in the parameter names. I think this also ensures backwards compatibility and would be preferred from my perspective.

Thanks.

Tagging @TannazH for information.

[rudymcintyre]

Hi Rudy,

Given your comment regarding gaussian_plume.py about passing arrays which could be different units I feel it might be a bit confusing to upate all the inner plume functions with more potential parameters which specify if it's degrees or m/s, leaving it like it is now is also a bit strange I think as it feels like mix of notation. So, I think the implementation you had before where the parameter was named the same only the unit might be different was a better way after all. Would it be possible to revert back to that commit, hence not have meter_per_sec or deg in the parameter names. I think this also ensures backwards compatibility and would be preferred from my perspective.

Thanks.

Tagging @TannazH for information.

Hey Bas,

Thanks for the feedback, that sounds good to me. Would you prefer if I revert to that commit? Or strip out the unit-specific bits but keep the documentation and other changes I made?

[bvandekerkhof]

I think it's probably easiets to strip out the unit specific bits. I really like what you did in the latest commit (Draxler Turbulence model, cbf5c54) with the turbulence model class, but probably the edits from the commit before that (Split turbulence vectors..., 6b3fe6c) need to be reverted. However going back to the commit before that also undo the latest edits which I would suggest to keep. Long story short: probably stripping the unit specific bits is the easiets. However, I'll leave it up to you what you prefer!