function [ n, e, u ] = Earth_deformation( fi0, la0, green_fun, Hydrology_model ) % % function Earth_deformation is used to calculate crustal deformations % as a response to hydrosphere loading % % fi0, la0 latitude and longitude of calculated deformations, % input as decimal degrees % % green_fun Green's function coefficients defining load response % of Earth's crust. % The data structure: 3 columns containing spherical distances % to loading mass and Green's function coefficients describing % deformations in radial and tangential direction. % Values in 2nd and 3rd columns should contain normalized % values (their de-normalization, depending on spherical % earth's radius is done at runtime). % Green's function coefficients for distances between the % distances in the first column of the table are calculated % with third degree spline % % Hydrology_model hydrological data with spatial resolution of % 0.5 degree (one month includes 360 rows and 720 columns) % that assign the amount of fresh water for each cell. % Data should express height of water column for each cell % in millimetres. Cells without determined water quantity % should be marked as NaN (Not a Number) %% additional variable format long; % load Green's function coefficients grn = green_fun; % load hydrology data to memory Hydro_model = Hydrology_model; % fresh water density ro = 1000; % kg/m^3 % name of file to store calculation results filename = 'Earth_def_fi_la_neu.txt'; %% parameters of reference surface % GRS80/WGS84 % a = 6378137; % e2 = 0.00669438002290; % spherical Earth a = 6371000; e2 = 0.0; %% Algorithm % fresh water column height at the point of calculated deformations H0 = hydro(Hydro_model, fi0, la0); % function calculates deformations only for land area and is skipping % any area that is missing WGHM data if isnan(H0) == 1 fprintf('Offshore location or WGHM data missing.\n') % if creating grid data for whole earth, % commend "end" below and % uncomment code between exclamation marks end %---!!!--- %{ n = 0; e = 0; u = 0; fileID = fopen(filename, 'at'); fprintf(fileID, '%2.1f %3.1f %2.8f %2.8f %2.8f\n', fi0, la0, n, e, u); fclose(fileID); else %} %---!!!--- dAz = 3; % degree(s) - width of spherical figures ds = 500; % basic integration distance [metres] R = a; r = ((pi*R)); % maximum integration distance - antipode of fi0, la0 nds = floor(r/ds); % amount of integration steps % coordinates of fi0, la0 antipode and fresh water column height a that location [fi5, la5] = directVincenty(fi0, la0, 0, pi*R, a, e2); H5 = hydro(Hydro_model, fi5, la5); % creation of empty tables in memory to store calculation results - improves execution time fig_area = zeros(nds+1,(360/dAz)); fig_hydro = zeros(nds+1,(360/dAz)); fig_mdistance = zeros(nds+1,(360/dAz)); fig_mazymut = zeros(nds+1, (360/dAz)); s = 1:nds; for Az = 0:(360/dAz)-1; Azymut = Az*dAz; %azimuth (bisector line) of spherical figure Az1 = Azymut-(dAz/2); Az2 = Azymut+(dAz/2); fig_mazymut(:,Az+1) = Azymut; % influence of one spherical triangle in neighbourhood of fi0, la0 s1 = s(1)*ds; [Area, H_hydro, dist] = triangle( fi0, la0, Az1, Az2, s1, a, e2, Hydro_model); fig_area(1,Az+1) = Area; fig_hydro(1,Az+1) = H_hydro; fig_mdistance(1,Az+1) = dist; % influence of one spherical trapezoid in zone s2 = (s(1,1:nds-1)*ds); s3 = (s(1,2:nds)*ds); [ Area, H_hydro, dist ] = trapezoid(fi0, la0, Az1, Az2, s2, s3, a, e2, Hydro_model); fig_area(2:end-1,Az+1) = Area'; fig_hydro(2:end-1,Az+1) = H_hydro'; fig_mdistance(2:end-1,Az+1) = dist'; % influence of one spherical triangle in neighbourhood of fi0, la0 antipode s4 = r - nds*ds; [Area, H_hydro, dist] = triangle( fi5, la5, Az1, Az2, s4, a, e2, Hydro_model); dist = pi*R - dist; fig_area(end,:) = Area; fig_hydro(end,:) = H_hydro; fig_mdistance(end,:) = dist; end % loading calculation, depending on spherical distance s = (grn(:,1)); L1_radial(:,1) = grn(:,2); L2_tangent(:,1) = grn(:,3); theta = deg2rad(s); %de-normalization of loading L1_radial = L1_radial ./ (a.*theta.*(10^12)); L2_tangent = L2_tangent ./ (a.*theta.*(10^12)); dist = fig_mdistance.*(180/(pi*a)); % interpolation of Green's function coefficients to match mean spherical distance % of spherical segment - 3rd degree spline G_radial = interp1(s,L1_radial,dist,'spline'); G_tangent = interp1(s,L2_tangent,dist,'spline'); % sum of individual segments on vertical deformations D_radial = G_radial.*fig_hydro.*fig_area; nany = isnan(D_radial); D_radial(nany) = 0; u = ro*sum(sum(D_radial)); % sum of individual segments on deformations % in direction of meridian and first vertical D_tangent = G_tangent.*fig_hydro.*fig_area; nany = isnan(D_tangent); D_tangent(nany) = 0; N = D_tangent.*(-cos(deg2rad(fig_mazymut))); E = D_tangent.*(-sin(deg2rad(fig_mazymut))); n = ro*sum(sum(N)); e = ro*sum(sum(E)); % saving calculated deformations to file fileID = fopen(filename, 'at'); fprintf(fileID, '%2.1f %3.1f %2.8f %2.8f %2.8f\n', fi0, la0, n, e, u); fclose(fileID); end

function [ N ] = Np( fi, a, e2) % % Np function is used to calculate radius of first vertical % at the latitude fi of given reference surface % % fi latitude for calculation of first vertical radius, % input as decimal degrees % % a radius of reference surface, input in metres % % e2 eccentricity of reference surface % % N radius of first vertical, % output as metres N = a./((1-e2.*((sin(deg2rad(fi))).^2)).^(1/2)); end

function [X, Y, Z] = blh2xyz(fi, la, h, a, e2) % % function blh2xyz calculate geocentric coordinates based on geodetic % coordinates and parameters of reference ellipsoid % % fi point's latitude, input as decimal degrees % % la point's longitude, input as decimal degrees % % a semi major axis of reference ellipsoid, input in metres % % e2 eccentricity of reference ellipsoid % % X, Y, Z geocentric coordinates of input point, output in metres %% Algorithm fi = deg2rad(fi); la = deg2rad(la); N = Np(fi, a, e2); X = ((N+h).*cos(fi).*cos(la)); Y = ((N+h).*cos(fi).*sin(la)); Z = (((N.*(1-e2))+h).*sin(fi)); end

function [y] = deg2rad(x) % % Function deg2rad converts degree angle to radians % % x input angle, enter as decimal degrees % % y output angle, expressed in radians y = x * pi / 180; end

function [ fi2, la2, a12, a21 ] = directVincenty( fi1, la1, alpha12, s12, a, e2) % % directVincenty function calculates coordinates of the geodesics end point % (fi2, la2) based on start point (fi1, la1), geodesic distance (s12) and % initial azimuth (alpha12) % % fi1 start point's latitude, input as decimal degrees % la1 start point's longitude, input as decimal degrees % alpha12 initial azimuth from start point to end point % (determined in start point), input as decimal degrees % s12 geodesic length, input as metres % % fi2 end point's latitude, output as decimal degrees % la2 end point's longitude, output as decimal degrees % alpha12 azimuth from start to end point, determined at the location % of end point, output as decimal degrees % alpha21 azimuth from end point to start point, determined at the location % of end point, output as decimal degrees % % Presented algorithm is based on direct Vincenty's formulae, to calculate % geodesics end point for reference surface and was altered to % improve calculation speed for spherical Earth %% convert decimal degrees input to radians fi1 = deg2rad(fi1); la1 = deg2rad(la1); a12 = deg2rad(alpha12); %% direct Vinceny's formulae % azimuth of geodesics intersecting the equator alpha = asin(cos(fi1)*sin(a12)); sigma = s12./a; fi2 = rad2deg(atan((sin(fi1).*cos(sigma)+cos(fi1).*sin(sigma).*cos(a12))./(((((sin(alpha)).^2)+(sin(fi1).*sin(sigma)-cos(fi1).*cos(sigma).*cos(a12)).^2).^(1/2))))); la = atan2(sin(sigma).*sin(a12),(cos(fi1).*cos(sigma)-sin(fi1).*sin(sigma).*cos(a12))); la2 = rad2deg(la1 + la); a12 = rad2deg(atan2(sin(alpha),(-sin(fi1).*sin(sigma)+cos(fi1).*cos(sigma).*cos(a12)))); if (a12 > -180 & a12 < 180) a21 = a12 + 180; else a21 = a12 - 180; end %fprintf('End point's coordinates:\nfi2 = %3.8f\nla2 = %3.8f\nazimuth12 = %3.8f\nazimuth21 = %3.8f\n',fi2, la2, a12, a21) end

Credits to Matsumoto K. “gotic2.”

0.0001 -33.64 -11.25 0.0010 -33.56 -11.25 0.0100 -32.75 -11.24 0.0200 -31.86 -11.21 0.0300 -30.98 -11.16 0.0400 -30.12 -11.09 0.0600 -28.44 -10.90 0.0800 -26.87 -10.65 0.1000 -25.41 -10.36 0.1600 -21.80 -9.368 0.2000 -20.02 -8.723 0.2500 -18.36 -8.024 0.3000 -17.18 -7.467 0.4000 -15.71 -6.725 0.5000 -14.91 -6.333 0.6000 -14.41 -6.150 0.8000 -13.69 -6.050 1.0 -13.01 -5.997 1.2 -12.31 -5.881 1.6 -10.95 -5.475 2.0 -9.757 -4.981 2.5 -8.519 -4.388 3.0 -7.533 -3.868 4.0 -6.131 -3.068 5.0 -5.237 -2.523 6.0 -4.660 -2.156 7.0 -4.272 -1.915 8.0 -3.999 -1.754 9.0 -3.798 -1.649 10.0 -3.640 -1.582 12.0 -3.392 -1.504 16.0 -2.999 -1.435 20.0 -2.619 -1.386 25.0 -2.103 -1.312 30.0 -1.530 -1.211 40.0 -0.292 -0.926 50.0 0.848 -0.592 60.0 1.676 -0.326 70.0 2.083 -0.223 80.0 2.057 -0.310 90.0 1.643 -0.555 100.0 0.920 -0.894 110.0 -0.025 -1.247 120.0 -1.112 -1.537 130.0 -2.261 -1.706 140.0 -3.405 -1.713 150.0 -4.476 -1.540 160.0 -5.414 -1.182 170.0 -6.161 -0.657 180.0 -6.663 0.000

function [ H ] = hydro( Hydro_model, fi, la ) % % function hydro is used to calculate water column height % at desired coordinates % % Hydro_model hydrological data with spatial resolution of 0.5 degree % (one month consists of 360 rows and 720 columns), with % information of water quantity as equivalent of fresh % water column with basis of 0.5 x 0.5 degrees and height % expressed in millimetres. % Data is expressed in millimetres. % No data is marked as NaN (Not a Number). % % fi, la coordinates for water column height determination, % input as decimal degrees % % water column height is determined by average sum of information % stored in four adjacent hydrosphere data cells %% Algorithm % Series of logical tests that returns numbers of columns and rows adjacent % to coordinates fi and la [vectorized] Urow = ceil(fi/0.5)+180'; Lrow = Urow; Rcol = ceil(la/0.5)'; zero = (Rcol == 0)*720; znak = sign(Rcol); znak(znak == 1) = 0; znak = znak.*-720; Rcol = Rcol+zero; Rcol = Rcol+znak; Lcol = Rcol; grid_fi = mod(fi,0.5) == 0; Urow = Urow + grid_fi; grid_la = mod(la,0.5) == 0; Rcol = Rcol + grid_la'; Urow = Urow'; Lrow = Lrow'; % returns index of elements adjacent to calculated deformations H_NE_idx = sub2ind(size(Hydro_model), Urow, Rcol); H_SE_idx = sub2ind(size(Hydro_model), Lrow, Rcol); H_SW_idx = sub2ind(size(Hydro_model), Lrow, Lcol); H_NW_idx = sub2ind(size(Hydro_model), Urow, Lcol); % returns water column height for every adjacent element H_NE = Hydro_model(H_NE_idx); H_SE = Hydro_model(H_SE_idx); H_SW = Hydro_model(H_SW_idx); H_NW = Hydro_model(H_NW_idx); [x] = [H_NE, H_SE, H_SW, H_NW]'; % calculate average height of water column, % sum and division skip Not a Number (NaN) %[H] = nansum(x)./(4-sum(isnan(x))); [H] = nanmean(x); end

function [ s12, alpha12, alpha21 ] = inverseVincenty( fi1, la1, fi2, la2, a, e2 ) % % inverseVincenty function is used to calculate geodesics distance s12 between two % points (fi1, fi1 and fi2, la2) on given reference surface (a, e2) as well % as azimuth of this geodesic % % fi1 start point's latitude, input as decimal degrees % la1 start point's longitude, input as decimal degrees % fi2 end point's latitude, input as decimal degrees % la2 end point's longitude, input as decimal degrees % a reference surface radius, input as metres % e2 reference surface eccentricity % % s12 geodesic length, output as metres % alpha12 azimuth from start point to end point, determined at the location % of start point, output as decimal degrees % alpha21 azimuth from end point to start point, determined at the location % of end point, output as decimal degrees % % %% additional calculations b = a*sqrt(1-e2); f = 1-sqrt(1-e2); fi1 = deg2rad(fi1); la1 = deg2rad(la1); fi2 = deg2rad(fi2); la2 = deg2rad(la2); if fi1 == fi2 && la1 == la2 fprintf('Start point and end point are the same\n') else %% inverse Vincenty's formulae % latitude reduced to sphere U1 = atan((1-f)*tan(fi1)); U2 = atan((1-f)*tan(fi2)); la = la2 - la1; L = la; for k=1:5 %k %L ssigma = sqrt((cos(U2).*sin(la)).^2+(cos(U1).*sin(U2)-sin(U1).*cos(U2).*cos(la)).^2); csigma = sin(U1).*sin(U2)+cos(U1).*cos(U2).*cos(la); sigma = atan2(ssigma,csigma); alpha = asin(cos(U1).*cos(U2).*sin(la)./ssigma); dwasigmam = acos(csigma-2.*sin(U1).*sin(U2)./((cos(alpha)).^2)); C = f./16.*((cos(alpha)).^2).*(4+f.*(4-3.*((cos(alpha)).^2))); la = L+(1-C).*f.*sin(alpha).*(sigma+C.*sin(sigma).*(cos(dwasigmam)+C.*cos(sigma).*(-1+2.*((cos(dwasigmam)).^2)))); %fprintf('dL = %1.16f \n', L-la); end u2 = ((cos(alpha)).^2).*(a.^2-b.^2)./(b.^2); A = 1+(u2./16384).*(4096+u2.*(-768+u2.*(320-175.*u2))); B = (u2./1024).*(256+u2.*(-128+u2.*(74-47.*u2))); % 1976b modification %k1 = (sqrt(1+u2)-1)/(sqrt(1+u2)+1); %A = (1+(k1^2)/4)/(1-k1); %B = k1*1-(3/8)*(k1^2); dsigma = B.*sin(sigma).*(cos(dwasigmam)+B./4.*(cos(sigma).*(-1+2.*((cos(dwasigmam)).^2)-B./6.*cos(dwasigmam).*(-3+4.*((sin(sigma)).^2)).*(-3+4.*((cos(dwasigmam)).^2))))); s12 = b.*A.*(sigma - dsigma); alpha12 = rad2deg(atan2(cos(U2).*sin(la),(cos(U1).*sin(U2)-sin(U1).*cos(U2).*cos(la)))); alpha21 = rad2deg(atan2(cos(U1).*sin(la),(-sin(U1).*cos(U2)+cos(U1).*sin(U2).*cos(la)))); %fprintf('Geodesic length s12 = %6.3f\nAzimuth12 = %3.8f \nAzimuth21 = %3.8f\n',s12, alpha12, alpha21); end

% script months is used to calculate deformations from hydrosphere % loading for desired coordinates for all hydrosphere data % approach uses spherical Earth as reference surface % fi0, la0 coordinates to calculate deformations, % input as decimal degrees % % green_function.mat Green's function coefficients, three columns: % 1st - spherical distances in decimal degrees % 2nd - coefficients for vertical direction % 3rd - coefficients for tangent direction % % WGHM.mat hydrosphere data with spatial resolution % of 0.5 x 0.5 degrees, expressed as height % (in millimetres) of water column in each cell % % Calculated deformation are expressed in millimetres and stored % in Earth_def_fi_la_neu.txt file. %warning ("off", "Octave:possible-matlab-short-circuit-operator") % Uncomment below line if Octave don't flush output to to console more off % load statistics package into Octave pkg load statistics clear delete Earth_def_fi_la_neu.txt % load Green's function coefficients and WGHM data to memory load grn1.mat load WGHM.mat % coordinates to calculate deformations fi0 = 52.1; la0 = 21.0; % loop divide hydrosphere data stored in WGHM.mat file to use in calculation for k = 1:length(WGHM)/360; fprintf('Month: %3i/%3i\n', k, length(WGHM)/360) month_model = WGHM(k*360-359:k*360,:); [ n, e, u ] = Earth_deformation(fi0, la0, grn1, month_model); end

function [ y ] = rad2deg( x ) % % rad2deg function converts radian angle to decimal degrees % % x input angle, input as radians % % y output angle, expressed in decimal degrees y = x / pi * 180; end

function [ Area, H_hydro, dist ] = trapezoid( fi0, la0, Az1, Az2, s1, s2, a, e2, Hydro_model ) % % trapezoid function is used to calculate properties of one spherical trapezoid % that influences hydrosphere loading % % fi0, la0 latitude and longitude of point for which deformation are calculated, % input as decimal degrees % % Az1, Az2 azimuths to left and right edges of spherical trapezoid, % input as decimal degrees % % s1, s2 distances to closer and further trapezoid bases, % (integration distances), input as metres % % a radius of reference surface, % input as metres % % e2 eccentricity of reference surface % % Hydro_model one epoch of hydrological data % % Area area of spherical trapezoid, expressed as square metres % % H_hydro average height of water column at spherical trapezoid, % expressed as millimetres % % dist distance to centre of spherical figure % %% calculation of figure's area % width of spherical trapezoid A = Az2-Az1; % calculate coordinates of spherical trapezoid's vertices [fi1, la1] = directVincenty(fi0, la0, Az1, s1, a, e2); [fi2, la2] = directVincenty(fi0, la0, Az2, s1, a, e2); [fi3, la3] = directVincenty(fi0, la0, Az1, s2, a, e2); [fi4, la4] = directVincenty(fi0, la0, Az2, s2, a, e2); % calculate spherical zone (spherical segment) height between % parallel planes, containing spherical trapezoid vertices, % perpendicular to sphere's rotation axis [fi5, la5] = directVincenty(90, 0, 0, s1, a, e2); [fi6, la6] = directVincenty(90, 0, 0, s2, a, e2); [~, ~, Z1] = blh2xyz(fi5, la5, 0, a, e2); [~, ~, Z2] = blh2xyz(fi6, la6, 0, a, e2); h = abs(Z2-Z1); Area = 2*pi*a*h*A/360; %% calculate average water column in spherical figure [H1] = hydro(Hydro_model,fi1,la1); [H2] = hydro(Hydro_model,fi2,la2); [H3] = hydro(Hydro_model,fi3,la3); [H4] = hydro(Hydro_model,fi4,la4); [x] = [H1; H2; H3; H4]; % calculate average water column, skipping Not-a-Number (NaN) H_hydro = nanmean(x); %% distance to centre of spherical figure dist = (s1+s2)/2; end

function [ Area, H_hydro, dist ] = triangle( fi0, la0, Az1, Az2, s, a, e2, Hydro_model ) % % triangle function is used to calculate properties of one spherical triangle % that influences hydrosphere loading % % fi0, la0 latitude and longitude of point for which deformation are calculated, % input as decimal degrees % % Az1, Az2 azimuths to left and right edges of spherical triangle, % input as decimal degrees % % s distance to vertices at Azimuth 1 and Azimuth 2 (integration distances), % input as metres % % a radius of reference surface, % input as metres % % e2 eccentricity of reference surface % % Hydro_model one epoch of hydrological data % % Area area of spherical triangle, expressed as square metres % % H_hydro average height of water column at spherical triangle, % expressed as millimetres % % dist distance to centre of spherical figure %% calculation of figure's area % width of spherical triangle A = Az2-Az1; % calculate coordinates of spherical triangle's vertices [fi1, la1] = directVincenty(fi0, la0, Az1, s, a, e2); %[fi1, la1, a01, a10] [fi2, la2] = directVincenty(fi0, la0, Az2, s, a, e2); %[fi2, la2, a02, a20] % use of inverse Vincenty's formulae to calculate distance between aforementioned % vertices [s12] = inverseVincenty(fi1, la1, fi2, la2, a, e2); %[s12, a12, a21] s12 = round(s12*10000)/10000; %rounds distance to 0.000001m % converts metric distance to radians s01 = deg2rad(s*360/(2*pi*a)); s12 = deg2rad(s12*360/(2*pi*a)); %[rad] % calculate spherical excess - L'Huilier's Theorem d = (s12+2*s01)/2; E = (atan(sqrt(tan(d./2).*tan((d-s12)./2).*tan((d-s01)./2).*tan((d-s01)./2)))).*4; Area = E*a^2; %% calculate average water column in spherical figure [H0] = hydro(Hydro_model,fi0,la0).*ones(size(fi1)); [H1] = hydro(Hydro_model,fi1(:),la1(:)); [H2] = hydro(Hydro_model,fi2(:),la2(:)); [x] = [H0; H1; H2]; % calculate average water column, skipping Not-a-Number (NaN) H_hydro = nanmean(x); %% distance to centre of spherical figure dist = s/2; end

%load WGHM_1.txt; all (all (WGHM == -9999)) %save('-mat7-binary','WGHM.mat','WGHM');