diff --git a/CHANGES.rst b/CHANGES.rst index 18be9c810..6114225f1 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -8,6 +8,7 @@ Contributors to this version: Éric Dupuis (:user:`coxipi`), Trevor James Smith New features and enhancements ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +* Add ``wind_power_potential`` to estimate the potential for wind power production given wind speed at the turbine hub height and turbine specifications, along with ``wind_profile`` to estimate the wind speed at different heights based on wind speed at a reference height. (:issue:`1458`, :pull:`1471`) * `xclim` now has a dedicated console command for prefetching testing data from `xclim-testdata` with branch options (e.g.: `$ xclim prefetch_testing_data --branch some_development_branch`). This command can be used to download the testing data to a local cache, which can then be used to run the testing suite without internet access or in "offline" mode. For more information, see the contributing documentation section for `Updating Testing Data`. (:issue:`1468`, :pull:`1473`). * The testing suite now offers a means of running tests in "offline" mode (using `pytest-socket `_ to block external connections). This requires a local copy of `xclim-testdata` to be present in the user's home cache directory and for certain `pytest` options and markers to be set when invoked. For more information, see the contributing documentation section for `Running Tests in Offline Mode`. (:issue:`1468`, :pull:`1473`). * The `SKIP_NOTEBOOKS` flag to speed up docs builds is now documented. See the contributing documentation section `Get Started!` for details. (:issue:`1470`, :pull:`1476`). @@ -35,6 +36,7 @@ Internal changes * GitHub deployment workflows now employs use of deployment environments for workflow security and uses the `Trusted Publisher `_ feature to sign and publish the `xclim` wheel and source distributions. (:pull:`1469`). * Mastodon publishing now uses `chuhlomin/render-template `_ and a standard formatting markdown document to format Mastodon toots. (:pull:`1469`). + v0.45.0 (2023-09-05) -------------------- Contributors to this version: David Huard (:user:`huard`), Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Juliette Lavoie (:user:`juliettelavoie`), Gabriel Rondeau-Genesse (:user:`RondeauG`), Marco Braun (:user:`vindelico`), Éric Dupuis (:user:`coxipi`). diff --git a/docs/references.bib b/docs/references.bib index e269a50ff..1ccf861a5 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -2046,3 +2046,30 @@ @misc{usda_2012 author = {{USDA Agricultural Research Service}}, year = {2012}, } + +@article{chen_2020, +title = {Impacts of climate change on wind resources over North America based on NA-CORDEX}, +journal = {Renewable Energy}, +volume = {153}, +pages = {1428-1438}, +year = {2020}, +issn = {0960-1481}, +doi = {https://doi.org/10.1016/j.renene.2020.02.090}, +url = {https://www.sciencedirect.com/science/article/pii/S0960148120302858}, +author = {Liang Chen}, +keywords = {Climate change, Wind energy, NA-CORDEX, North America} +} + +@article{tobin_2018, +doi = {10.1088/1748-9326/aab211}, +url = {https://dx.doi.org/10.1088/1748-9326/aab211}, +year = {2018}, +month = {apr}, +publisher = {IOP Publishing}, +volume = {13}, +number = {4}, +pages = {044024}, +author = {I Tobin and W Greuell and S Jerez and F Ludwig and R Vautard and M T H van Vliet and F-M Bréon}, +title = {Vulnerabilities and resilience of European power generation to 1.5°C, 2°C and 3°C warming}, +journal = {Environmental Research Letters} +} diff --git a/tests/test_atmos.py b/tests/test_atmos.py index 2828f92fc..5c86b8c29 100644 --- a/tests/test_atmos.py +++ b/tests/test_atmos.py @@ -253,6 +253,41 @@ def test_wind_chill_index(atmosds): ) +def test_wind_profile(atmosds): + out = atmos.wind_profile( + wind_speed=atmosds.sfcWind, h_r="10 m", h="100 m", alpha=1 / 7 + ) + assert out.attrs["units"] == "m s-1" + assert (out > atmosds.sfcWind).all() + + +def test_wind_power_potential(atmosds): + out = atmos.wind_power_potential(wind_speed=atmosds.sfcWind) + assert out.attrs["units"] == "" + assert (out >= 0).all() + assert (out <= 1).all() + + +def test_wind_power_potential_from_3h_series(): + """Test a typical computation workflow from 3-hourly time series to daily power production in MWh.""" + from xclim.core.units import convert_units_to + from xclim.indices.generic import select_resample_op + from xclim.testing.helpers import test_timeseries + + w = test_timeseries( + np.ones(96) * 15, variable="sfcWind", start="7/1/2000", units="m s-1", freq="3H" + ) + out = atmos.wind_power_potential(wind_speed=w) + + # Multiply with nominal capacity + power = out * 100 + power.attrs["units"] = "MW" + annual_power = convert_units_to( + select_resample_op(power, op="sum", freq="D"), "MWh" + ) + assert (annual_power == 100 * 24).all() + + class TestDrynessIndex: def test_simple(self, atmosds): ds = atmosds.isel(location=3) diff --git a/tests/test_indices.py b/tests/test_indices.py index 3f7176a38..613669582 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -3626,3 +3626,29 @@ def test_late_frost_days_lat(self, tasmin_series): tasmin = tasmin_series(np.array([-1, 1, 2, -4, 0]) + K2C, start="30/3/2023") lfd = xci.frost_days(tasmin, date_bounds=("04-01", "06-30")) np.testing.assert_allclose(lfd, 1) + + +class TestWindProfile: + def test_simple(self, sfcWind_series): + a = np.linspace(0, 100) + v = xci.wind_profile(sfcWind_series(a), h="100 m", h_r="10 m") + np.testing.assert_allclose(v, a * 10 ** (1 / 7)) + + +class TestWindPowerPotential: + def test_simple(self, sfcWind_series): + v = [2, 6, 20, 30] + p = xci.wind_power_potential( + sfcWind_series(v, units="m/s"), cut_in="4 m/s", rated="8 m/s" + ) + np.testing.assert_allclose(p, [0, (6**3 - 4**3) / (8**3 - 4**3), 1, 0]) + + # Test discontinuities at the default thresholds + v = np.array([3.5, 15]) + a = sfcWind_series(v - 1e-7, units="m/s") + b = sfcWind_series(v + 1e-7, units="m/s") + + pa = xci.wind_power_potential(a) + pb = xci.wind_power_potential(b) + + np.testing.assert_array_almost_equal(pa, pb, decimal=6) diff --git a/xclim/core/units.py b/xclim/core/units.py index b79fddb08..474394afb 100644 --- a/xclim/core/units.py +++ b/xclim/core/units.py @@ -579,7 +579,7 @@ def to_agg_units( units="", is_dayofyear=np.int32(1), calendar=get_calendar(orig) ) - elif op in ["count", "integral"]: + elif op in ["count", "integral", "sum"]: m, freq_u_raw = infer_sampling_units(orig[dim]) orig_u = str2pint(orig.units) freq_u = str2pint(freq_u_raw) @@ -587,7 +587,7 @@ def to_agg_units( if op == "count": out.attrs["units"] = freq_u_raw - elif op == "integral": + elif op in ["integral", "sum"]: out.attrs["units"] = pint2cfunits(orig_u * freq_u) else: raise ValueError( diff --git a/xclim/data/fr.json b/xclim/data/fr.json index 646a1f6a9..96356ae89 100644 --- a/xclim/data/fr.json +++ b/xclim/data/fr.json @@ -792,6 +792,18 @@ "long_name": "Indice forêt météo", "description": "Évaluation numérique de l'intensité du feu." }, + "WIND_POWER_POTENTIAL": { + "long_name": "Potentiel de production éolienne", + "description": "Potentiel de production éolienne estimé par une courbe de puissance semi-idéalisée d'une turbine, paramétrée par une vitesse minimale de démarrage {cut_in}, la vitesse nominale {rated}, et la vitesse d'arrêt {cut_out}.", + "title": "Potentiel de production éolienne", + "abstract": "Estimation à partir d'une loi de puissance semi-idéalisée de la production d'énergie éolienne selon la vitesse du vent." + }, + "WIND_PROFILE": { + "long_name": "Vitesse du vent à la hauteur {h}", + "description": "Vitesse du vent à une hauteur de {h} calculée de la vitesse à {h_r} par un profil de loi de puissance.", + "title": "Profil du vent", + "abstract": "Estimation de la vitesse du vent à partir du vent à une hauteur de référence." + }, "WIND_SPEED_FROM_VECTOR": { "title": "Vitesse et direction du vent à partir du vecteur vent", "abstract": "Calcul de la magnitude et la direction de la vitesse du vent à partir des deux composantes ouest-est et sud-nord." diff --git a/xclim/data/variables.yml b/xclim/data/variables.yml index 4f14d5e3b..0a7553b0c 100644 --- a/xclim/data/variables.yml +++ b/xclim/data/variables.yml @@ -1,4 +1,10 @@ variables: + air_density: + canonical_units: kg m-3 + cell_methods: "time: mean" + description: Air density. + dimensions: "[density]" + standard_name: air_density areacello: canonical_units: m2 cell_methods: "area: sum" @@ -316,6 +322,12 @@ variables: description: Northward component of the wind velocity (at the surface). dimensions: "[speed]" standard_name: northward_wind + wind_speed: + canonical_units: m s-1 + cell_methods: "time: mean" + description: Wind speed. + dimensions: "[speed]" + standard_name: wind_speed wsgsmax: cmip6: False canonical_units: m s-1 diff --git a/xclim/indicators/atmos/_conversion.py b/xclim/indicators/atmos/_conversion.py index d1771e88f..6876755f2 100644 --- a/xclim/indicators/atmos/_conversion.py +++ b/xclim/indicators/atmos/_conversion.py @@ -28,6 +28,8 @@ "water_budget", "water_budget_from_tas", "wind_chill_index", + "wind_power_potential", + "wind_profile", "wind_speed_from_vector", "wind_vector_from_speed", ] @@ -123,6 +125,30 @@ def cfcheck(self, **das): compute=indices.sfcwind_2_uas_vas, ) +wind_power_potential = Converter( + title="Wind power potential", + identifier="wind_power_potential", + units="", + long_name="Wind power potential", + description="Wind power potential using a semi-idealized turbine power curve using a cut_in speed of {cut_in}, " + "a rated speed of {rated}, and a cut_out speed of {cut_out}.", + cell_methods="", + abstract="Calculation of the wind power potential using a semi-idealized turbine power curve.", + compute=indices.wind_power_potential, +) + +wind_profile = Converter( + title="Wind profile", + identifier="wind_profile", + var_name="wind_speed", + units="m s-1", + standard_name="wind_speed", + long_name="Wind speed at height {h}", + description="Wind speed at a height of {h} computed from the wind speed at {h_r} using a power law profile.", + cell_methods="", + abstract="Calculation of the wind speed at a given height from the wind speed at a reference height.", + compute=indices.wind_profile, +) saturation_vapor_pressure = Converter( title="Saturation vapour pressure (e_sat)", diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py index 5183c6531..6e571e082 100644 --- a/xclim/indices/_conversion.py +++ b/xclim/indices/_conversion.py @@ -53,6 +53,8 @@ "uas_vas_2_sfcwind", "universal_thermal_climate_index", "wind_chill_index", + "wind_power_potential", + "wind_profile", ] @@ -2017,3 +2019,161 @@ def mean_radiant_temperature( 0.25, ) return mrt.assign_attrs({"units": "K"}) + + +@declare_units(wind_speed="[speed]", h="[length]", h_r="[length]") +def wind_profile( + wind_speed: xr.DataArray, + h: Quantified, + h_r: Quantified, + method: str = "power_law", + **kwds, +): + r"""Wind speed at a given height estimated from the wind speed at a reference height. + + Estimate the wind speed based on a power law profile relating wind speed to height above the surface. + + Parameters + ---------- + wind_speed : xarray.DataArray + Wind speed at the reference height. + h : Quantified + Height at which to compute the wind speed. + h_r : Quantified + Reference height. + method : {"power_law"} + Method to use. Currently only "power_law" is implemented. + kwds : dict + Additional keyword arguments to pass to the method. For power_law, this is alpha, which takes a default value + of 1/7, but is highly variable based on topography, surface cover and atmospheric stability. + + Notes + ----- + The power law profile is given by + + .. math:: + + v = v_r \left( \frac{h}{h_r} \right)^{\alpha}, + + where :math:`v_r` is the wind speed at the reference height, :math:`h` is the height at which the wind speed is + desired, and :math:`h_r` is the reference height. + + """ + # Convert units to meters + h = convert_units_to(h, "m") + h_r = convert_units_to(h_r, "m") + + if method == "power_law": + alpha = kwds.pop("alpha", 1 / 7) + out = wind_speed * (h / h_r) ** alpha + out.attrs["units"] = wind_speed.attrs["units"] + return out + else: + raise NotImplementedError(f"Method {method} not implemented.") + + +@declare_units( + wind_speed="[speed]", + air_density="[air_density]", + cut_in="[speed]", + rated="[speed]", + cut_out="[speed]", +) +def wind_power_potential( + wind_speed: xr.DataArray, + air_density: xr.DataArray = None, + cut_in: Quantified = "3.5 m/s", + rated: Quantified = "13 m/s", + cut_out: Quantified = "25 m/s", +) -> xr.DataArray: + r"""Wind power potential estimated from an idealized wind power production factor. + + The actual power production of a wind farm can be estimated by multiplying its nominal (nameplate) capacity by the + wind power potential, which depends on wind speed at the hub height, the turbine specifications and air density. + + Parameters + ---------- + wind_speed : xarray.DataArray + Wind speed at the hub height. Use the `wind_profile` function to estimate from the surface wind speed. + air_density: xarray.DataArray + Air density at the hub height. Defaults to 1.225 kg/m³. This is worth changing if applying in cold or + mountainous regions with non-standard air density. + cut_in : Quantified + Cut-in wind speed. Default is 3.5 m/s. + rated : Quantified + Rated wind speed. Default is 13 m/s. + cut_out : Quantified + Cut-out wind speed. Default is 25 m/s. + + Returns + ------- + xr.DataArray + The power production factor. Multiply by the nominal capacity to get the actual power production. + + See Also + -------- + wind_profile : Estimate wind speed at the hub height from the surface wind speed. + + Notes + ----- + This estimate of wind power production is based on an idealized power curve with four wind regimes specified + by the cut-in wind speed (:math:`u_i`), the rated speed (:math:`u_r`) and the cut-out speed (:math:`u_o`). + Power production is zero for wind speeds below the cut-in speed, increases cubically between the cut-in + and rated speed, is constant between the rated and cut-out speed, and is zero for wind speeds above the cut-out + speed to avoid damage to the turbine :cite:p:`tobin_2018`: + + .. math:: + + \begin{cases} + 0, & v < u_i \\ + (v^3 - u_i^3) / (u_r^3 - u_i^3), & u_i ≤ v < u_r \\ + 1, & u_r ≤ v < u_o \\ + 0, & v ≥ u_o + \end{cases} + + For non-standard air density (:math:`\rho`), the wind speed is scaled using + :math:`v_n = v \left( \frac{\rho}{\rho_0} \right)^{1/3}`. + + The temporal resolution of wind time series has a significant influence on the results: mean daily wind + speeds yield lower values than hourly wind speeds. Note however that percent changes in the wind power potential + climate projections are similar across resolutions :cite:p:`chen_2020`. + + To compute the power production, multiply the power production factor by the nominal + turbine capacity (e.g. 100), set the units attribute (e.g. "MW"), resample and sum with + `xclim.indices.generic.select_resample_op(power, op="sum", freq="D")`, then convert to + the desired units (e.g. "MWh") using `xclim.core.units.convert_units_to`. + + References + ---------- + :cite:cts:`chen_2020,tobin_2018`. + + """ + # Convert units + cut_in = convert_units_to(cut_in, wind_speed) + rated = convert_units_to(rated, wind_speed) + cut_out = convert_units_to(cut_out, wind_speed) + + # Correct wind speed for air density + if air_density is not None: + default_air_density = convert_units_to("1.225 kg/m^3", air_density) + f = (air_density / default_air_density) ** (1 / 3) + else: + f = 1 + + v = wind_speed * f + + out = xr.apply_ufunc(_wind_power_factor, v, cut_in, rated, cut_out) + out.attrs["units"] = "" + return out + + +@vectorize +def _wind_power_factor(v, cut_in, rated, cut_out): + """Wind power factor function""" + if v < cut_in: + return 0.0 + if v < rated: + return (v**3 - cut_in**3) / (rated**3 - cut_in**3) + if v < cut_out: + return 1.0 + return 0.0