pvl_soiling_hsu.m

function SR = pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, varargin)
% PVL_SOILING_HSU Calculates soiling rate over time given particulate and rain data
%
% Syntax
%   [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10)
%   [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, ModelType)
%   [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, ModelType, RainAccPeriod)
%   [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, ModelType, RainAccPeriod, LUC, WindSpeed, Temperature)
%
% Description  
%   [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10)
%      Uses the time, rainfall amounts, tilt angle, particulate matter
%      concentrations, and rain cleaning threshold to predict soiling rate
%      using the fixed settling model. Uses hourly rain accumulation to
%      determine cleaning.
%   [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, ModelType)
%      Uses the time, rainfall amounts, tilt angle, particulate matter
%      concentrations, and rain cleaning threshold to predict soiling rate
%      using the model type selected by ModelType. Uses hourly rain
%      accumulation to determine cleaning.
%   [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, ModelType, RainAccPeriod)
%      Uses the time, rainfall amounts, tilt angle, particulate matter
%      concentrations, and rain cleaning threshold to predict soiling rate
%      using the model type selected by ModelType. The user specifies the
%      rain accumulation period (in hours) over which to determine whether
%      sufficient rain has fallen to clean the module.
%   [SR]=pvl_soiling_hsu(Time, Rain, RainThresh, Tilt, PM2_5, PM10, ModelType, RainAccPeriod, LUC, WindSpeed, Temperature)
%      Uses the parameters of the fixed models along with Land Use Category
%      (LUC), wind speed, and ambient temperature to predict soiling ratio
%      using the variable deposition velocity model.
%
% Input Parameters:
%   Time is a struct with elements as described in pvl_maketimestruct. 
%     Time values for the soiling function do not need to be 
%     regularly spaced, although large gaps in timing are discouraged.
%
%   Rain is a vector of rainfall in mm of the same length as Time, where
%     each element of Rain correlates with the element of Time.
%     Rainfall values should be in mm of rainfall. Programmatically, rain
%     is accumulated over a given time period, and cleaning is applied
%     immediately after a time period where the cleaning threshold is
%     reached.
%
%   RainThresh is a scalar for the amount of rain, in mm, in an accumulation
%     period needed to clean the PV modules. In periods where the 
%     accumulated rain meets or exceeds RainThresh, the panels are assumed 
%     to be cleaned immediately after the accumulation period 
%     [1] suggests a good RainThresh could be 1mm, but the time period is
%     not specified. Rain accumulation period length can be adjusted in the
%     optional input RainAccPeriod.
%
%   Tilt is a scalar or vector for the tilt of the PV panels. Changing tilt
%     angles (e.g. in tracking cases) can be accomodated, and tilt angles
%     are correlated with the entries in Time.
%
%   PM2_5 is the concentration of airborne particulate matter (PM) with diameter less 
%     than 2.5 microns, in g/m^3. Historical US concentration data sets may be
%     found at https://aws.epa.gov/aqsweb/airdata/download_files.html. With
%     file descriptions (for hourly data sets) at 
%     https://aqs.epa.gov/aqsweb/airdata/FileFormats.html#_hourly_data_files
%   
%   PM10 is the concentration of airborne particulate matter (PM) with diameter less than
%     10 microns, in g/m^3. Historical US PM data sets found as described
%     above.
%
%   ModelType is an optional input to the function to determine the 
%     the model type to be used in the soiling model, see [1]. A
%     value of "1" indicates that the Variable Deposition Velocity model
%     shall be used, a value of "2" indicates that the Fixed Settling
%     Velocity model shall be used, and a value of "3" indicates that the
%     Fixed Deposition Velocity model shall be used. [1] indicates that the
%     Fixed Settling Velocity model performs best under a wide range of
%     conditions, and thus "2" is the default ModelType if ModelType is omitted. 
%     Validation efforts by Sandia National Laboratories
%     confirm these findings. If an incorrect ModelType is provided, the
%     Fixed Settling Velocity (type 2) will be used (with a warning).
%
%   RainAccPeriod is an optional input that specifies the period, in hours,
%     over which to accumulate rainfall totals before checking against the
%     rain cleaning threshold. For example, if the rain threshold is
%     0.5 mm per hour, then RainThresh should be 0.5 and RainAccPeriod
%     should be 1. If the threshold is 1 mm per hour, then the values
%     should be 1 and 1, respectively. The minimum RainAccPeriod is 1
%     hour. The default value is 1, indicating hourly rain accumulation.
%     Accumulation periods exceeding 24 (daily accumulation) are not
%     recommended.
%   
%   LUC is an optional input to the function, but it is required for the
%     Variable Deposition Model. LUC is the Land Use Category as specified
%     in Table 19.2 of [2]. LUC must be a numeric scalar with value 1, 4,
%     6, 8, or 10, corresponding to land with evergreen trees, deciduous
%     trees, grass, desert, or shrubs with interrupted woodlands. If 
%     omitted, the default value of 8 (desert) is used. 
%
%   WindSpeed is an optional input to the function, but is required for the
%     Variable Deposition Model. WindSpeed is a scalar or vector value with
%     the same number of elements as Time, and must be in meters per
%     second. If WindSpeed is omitted, the value of 2 m/s is used as
%     default.
%
%   Temperature is an optional input to the function, but is required for
%     the Variable Deposition Model. Temperature is a scalar or vector
%     value with the same number of Elements as Time and must be in degrees
%     C. If Temperature is omitted, the value of 12 C is used as default.
% 
% 
% Output Parameters:
%   SR = The soiling ratio (SR) of a tilted PV panel, this is a number
%   between 0 and 1. SR is a time series where each element of SR
%   correlates with the accumulated soiling and rain cleaning at the times
%   specified in Time.
%
% Background and recognition:
%   This soiling model was developed by Liza Boyle and Merissa Coello at
%   Humboldt State University (HSU). Many thanks to Liza and Merissa for
%   model development, initial code, and subsequent code and performance
%   reviews.
%
% References
%   [1] M. Coello and L. Boyle, "Simple Model For Predicting Time Series 
%       Soiling of Photovoltaic Panels," in IEEE Journal of Photovoltaics.
%       doi: 10.1109/JPHOTOV.2019.2919628
%   [2] Atmospheric Chemistry and Physics: From Air Pollution to Climate
%       Change. J. Seinfeld and S. Pandis. Wiley and Sons 2001.
%
% See also PVL_MAKETIMESTRUCT 

p=inputParser;

p.addRequired('Time', @isstruct);
p.addRequired('Rain', @(x) all(isnumeric(x) & all(x>=0) & isvector(x)));
p.addRequired('RainThresh', @(x) all(isnumeric(x) & isscalar(x) & all(x>=0)));
p.addRequired('Tilt', @(x) all(isnumeric(x) & isvector(x)));
p.addRequired('PM2_5', @(x) all(isnumeric(x) & all(x>=0) & isvector(x)))
p.addRequired('PM10', @(x) all(isnumeric(x) & all(x>=0) & isvector(x)))
p.addOptional('ModelType', 2, @(x) all(isscalar(x)));
p.addOptional('RainAccPeriod', 1, @(x) all(isscalar(x) & x>=1 & isnumeric(x)));
p.addOptional('LUC', 8, @(x) all(isscalar(x)))
p.addOptional('WindSpeed', 2, @(x) all(isnumeric(x) & all(x>=0) & isvector(x)));
p.addOptional('Temperature', 12, @(x) all(isnumeric(x) & all(x>=0) & isvector(x)));
p.parse(Time, Rain, RainThresh, Tilt, PM2_5, PM10, varargin{:});

RainAccPeriod = p.Results.RainAccPeriod;
ModelType = p.Results.ModelType;
LUC = p.Results.LUC;
WindSpeed = p.Results.WindSpeed;
Temperature = p.Results.Temperature;

% Create a datenum from the time structure to handle duration easily
TimeAsDatenum = datenum(Time.year, Time.month, Time.day, Time.hour, ...
    Time.minute, Time.second);

% Following section accumulates the rain in each accumulation period
% To change this to a more frequent accumulation rate, create a
% new array based on TimeAsDatenum where the unit is the given unit of
% accumulation (e.g. for hourly accumulation you might use
% TimeAsDatenum*24)
RainAccAsDatenum = floor(TimeAsDatenum * 24 ./ RainAccPeriod);
% Determine the number of unique days, the first instance of that day, and
% create an index of which entries are part of that day
[RainAccTimes, UnqRainAccFrstVal, UnqRainAccIndx] = unique(RainAccAsDatenum);
% Accumulate rain readings according to the unique index of days
RainAtAccTimes = accumarray(UnqRainAccIndx, Rain);
% Create an accumulated rain vector, then populate the last entry of each
% date with the rain that fell on that day
AccumRain = zeros(size(Rain));
AccumRain(UnqRainAccFrstVal(2:end)-1) = RainAtAccTimes(1:end-1);
% Rain accumulated on the last day is placed at the last time (does not
% affect soiling calculation)
AccumRain(end) = RainAtAccTimes(end);

switch ModelType
    case 1 % Variable Deposition Velocity
        vd = depo_veloc(Temperature, WindSpeed, LUC);
    case 2 % Fixed Settling Velocity in m/s
        vd(1,2) = 0.004;
        vd(1,1) = 0.0009;
    case 3 % Fixed Deposition Velcoity in m/s
        vd(1,2) = 0.0917;
        vd(1,1) = 0.0015;
    otherwise
        warning('Unknown Model Type specified, using Fixed Setting Velocity model')
        vd(1,2) = 0.004;
        vd(1,1) = 0.0009;
end


% Corrects PM10 measurement since PM2.5 is included in measurement
PMConcentration=nan(numel(TimeAsDatenum),2); % pre-allocate with NAN

PMConcentration(:,1) = PM2_5; % fill PM2.5 data in column 1
% Since PM2.5 particulate is counted in PM10 measurements, and we need only
% the particular from 2.5 micron to 10 micron, subtract the 2.5 micron
% particulate from the 10 micron particulate
PMConcentration(:,2) = PM10 - PM2_5; % fill in PM2.5-PM10 data in column 2

% For any instances where 2.5 micron particulate was more thatn 10 micron,
% set the amount of 2.5 to 10 micron particulate to 0
PMConcentration(PM10 - PM2_5 < 0 , 2) = 0; 

% EPA concentrations reported in micrograms/m^3, we need it in g/m^3
PMConcentration = PMConcentration .* 1e-6;

% Calculate the mass deposition rate of each particulate size
F=PMConcentration .* vd; % g * m^-2 * s^-1, by particulate size
HorizontalTotalMassRate = F(:,1) + F(:,2); % g * m^-2 * s^-1, total

% Mass depositition rate on a tilted surface. In the event of a changing
% tilt (e.g. tracking system) the mass rate is adjusted by the tilt angle
% at each time step, and thus a linear interpolation of changing tilt
% angles are used (when using trapezoidal integration).
TiltedMassRate = HorizontalTotalMassRate .* cosd(Tilt);

% Mass on the tilted surface if no rain occurs. Uses trapezoidal
% integration and determines the amount of mass deposited on a tilted
% module.
TiltedMassNoRain = cumtrapz(TimeAsDatenum*86400, TiltedMassRate);

% Create a Tilted soiling mass that will account for rain cleaning
TiltedMass = TiltedMassNoRain;
% For each rain accumulation period, check to see if enough rain fell to meet or
% exceed the cleaning threshold. If the threshold was met or exceeded,
% then subtract from each subsuequent time step the mass of future time
% steps. This essentially "resets" the cumulative soiling accumulation when
% the rain meets or exceeds the stated thresshold. Note that this means
% that the mass will never be reported as zero because the cleaning is
% applied immediately after a time step, and new mass accumulates between
% the cleaning and the first time reported after the cleaning.
% Note that this could be done more clearly by going through each time step
% and accumulating or clearing based on daily rainfall, but this is much
% faster, and can utilize trapezoidal soiling integration.
for cntr1 = 1:numel(RainAtAccTimes)-1
    if RainAtAccTimes(cntr1) >= RainThresh
        TiltedMass(UnqRainAccFrstVal(cntr1+1):end) = TiltedMass(UnqRainAccFrstVal(cntr1+1):end)-TiltedMass(UnqRainAccFrstVal(cntr1+1)-1);
    end
end

% Soiling rate is the % transmittance reduction using Heagzy equation
SoilingRate = 34.37 * erf(0.17 .* (TiltedMass.^0.8473));

% Soiling Ratio is the % soil transmittance (0.95 = 5% soiling loss)
SR = (100 - SoilingRate) ./ 100;

end

% This function creates deposition velocities for the variable deposition
% model.
function vd = depo_veloc(T, WindSpeed, LUC)

    % convert temperature into Kelvin 
    T = T + 273.15;

    % save wind data
    u = WindSpeed;

    g=9.81;         %gravity in m/s^2
    Na=6.022*10^23; %avagadros number
    R=8.314;        %Universal gas consant in m3Pa/Kmol
    k=1.38*10^-23;  %Boltzmann's constant in m^2kg/sK
    P=101300;       %pressure in Pa
    rhoair= 1.2041; %density of air in kg/m3
    z0=1;
    rhop=1500;      %Assume density of particle in kg/m^3
    
    switch LUC
        case {1, 4}
            gamma = 0.56;
        case {6, 8, 10}
            gamma = 0.54;
        otherwise
            warning('Unknown Land Use Category (LUC), assuming LUC 8.');
            gamma = 0.54;
    end


    %Diameter of particle in um
    Dpum=[2.5,10];
    Dpm=Dpum*10^-6;   %Diameter of particle in m

    %Calculations
    mu=1.8*10^-5.*(T./298).^0.85;      %viscosity of air in kg/m s
    nu=mu/rhoair;
    lambda1=2*mu./(P.*(8.*0.0288./(pi.*R.*T)).^(1/2));   %mean free path
    
    Cc = 1+2.*[lambda1 lambda1]./Dpm.*(1.257+0.4.*exp(-1.1.*Dpm./([lambda1 lambda1].*2))); %slip correction coefficient

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Calculate vs
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    vs = rhop.*Dpm.^2 .* (g .* Cc ./(mu.*18)); %particle settling velocity

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Calculate rb
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ustar = NaN(size(u)); % pre-allocate ustar
    %Equation 11.66 in Ramaswami (and 16.67 and Sienfeld &Pandis)
    ustar(u > 0) = 0.4 * u(u > 0) ./ log(10/z0);
    ustar(u<=0) = 0.001;

    D=k*T.*(Cc./(3*pi*mu*Dpm));

    Sc=nu./D;
    %gamma=.56;      %for urban
    %alpha=1.5;     %for urban      
    EB=Sc.^(-1 * gamma);
    St = vs .* (ustar .^ 2) ./ (g .* nu);
    
    EIM=10.^(-3./St);   %For smooth surfaces
    %EIM=((St)./(0.82+St)).^2;

    R1=exp(-St.^(1/2));  %percentage of particles that stick

    rb=1./(3*(EB+EIM).*ustar.*R1);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Calculate ra
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    a=[-0.096,-0.037,-0.002,0,0.004,0.035];
    b=[0.029,0.029,0.018,0,-0.018,-0.036];
    
    % For wind speeds <= 3, use a = -0.037 and b = 0.029
    % For wind speeds >3 and <=5, use a = -.002, b = 0.018
    % For wind speeds > 5, use a = 0, b = 0
    avals = a(2) .* ones(size(u));
    avals(u>3) = a(3);
    avals(u>5) = a(4);
    
    bvals = b(2) .* ones(size(u));
    bvals(u>3) = b(3);
    bvals(u>5) = b(4);
    
    L = 1 ./ (avals + bvals .* log(z0));
    
    zeta0=z0./L;
    zeta=10./L;
    eta = ((1-15.*zeta).^(1./4));
    eta0 = ((1-15.*zeta0).^(1./4));
    
    ra = NaN(size(zeta)); % Preallocate memory
    ra(zeta == 0) = (1 ./ (0.4 .* ustar(zeta == 0))) .* log(10 ./ z0);
    ra(zeta > 0) = (1./(0.4.*ustar(zeta > 0))).*(log(10./z0) + ...
        4.7.*(zeta(zeta > 0)-zeta0(zeta > 0)));
    ra(zeta < 0) = (1 ./ (0.4 .* ustar(zeta < 0))) .* (log(10 ./ z0) + ...
        log((eta0(zeta < 0).^2 + 1) .* (eta0(zeta < 0)+1).^2 ./ ...
        ((eta(zeta < 0).^2 + 1) .* (eta(zeta < 0)+1).^2)) + ...
        2 .* (atan(eta(zeta < 0))-atan(eta0(zeta < 0))));


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Calculate vd and mass flux
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    vd=1./(ra+rb)+vs;

end