Model Description
The authors address the question of how infrastructure design affects SES sustainability in two stages. First, they explore the effects of design variations in shared infrastructure on long-term system behavior in a model system. They examine two types of distribution infrastructure, one with and one without upstream−downstream asymmetry, and different threshold characteristics of infrastructure maintenance. Second, they evaluate how these design variations influence the robustness of system function to an economic shock.
The model assumes that there are N farming households spread across two villages (village 1 and village 2) that manage a shared irrigation infrastructure. Each farmer is endowed with the same amount of available labor (l) each year and the same acreage (a). A farmer appropriates q units of water from the system and allocates labor among three activities: farming, maintaining infrastructure and outside employment at a wage rate w. Governance is represented in the model by the following rules. The expected maintenance labor contribution is proportional to the farmer's acreage (assumed to be same for all farmers in this model). Water allocations are also proportional to acreage, but only among the water rights holders. Only farmers who contributed labor to the infrastructure maintenance before the planting season obtain water rights. Farmers choose between two strategies: group conformist (G) and opportunist (O). Gs follow and enforce rules, and strive to maximize the total welfare of the two villages. Os break the rules and attempt to maximize the individual net income. The model tracks the fraction of Gs in the village i and the resulting performance of the irrigation system.
$\begin{equation}
\frac{dX_i}{dt} = X_i(\pi_i^G - \overline{\pi_i})
\end{equation}$ |
| This replicator equation tracks the fractions of the strategies in villages 1 and 2. |
$\begin{equation}
Q = I(L_m)S(t)
\end{equation}$ |
| The total irrigation water, Q, depends on the efficiency of the shared irrigation infrastructure and the volume of river discharge. Efficiency of the irrigation infrastructure depends on the sum of maintenance labor provided by all farmers each year. |
$I(L_m) =
\begin{cases}
0 & 0 \leq L_m < \psi - \epsilon \\
\frac{I_{max}}{2\epsilon} (L_m - \psi + \epsilon) & \psi - \epsilon \leq L_m \leq \psi + \epsilon\\
I_{max} & \psi + \epsilon \leq L_m \leq L
\end{cases}$ |
| I-max is the maximum infrastructure efficiency. The half-saturation point of labor and half-width of the threshold slope determine the threshold of maintenance. |
$\begin{equation}
l = l_f + l_m + l_e
\end{equation}$ |
| A farmer can allocate labor among three activities: farming (l-f), maintaining infrastructure (l-m) and outside employment (l-e) |
$\begin{equation}
\pi_i = pb(l_f)^j(q_i + r_i)^ka_i^(1-j-k) + wl_e
\end{equation}$ |
| The equation represents a farmer's total income in village i |
$\begin{equation}
\begin{aligned}
& \underset{L_f,L_m}{\text{max}}
{\Pi = pbL_f^j[I(L_m)S + R]^k A^(1-j-k) + w L_e} \\
&\text{s.t.} L_f + L_m + L_e = L \\
\end{aligned}
\end{equation}$ |
| The aggregate income of two villages maximized over the optimal values of total maintenance labor, farming, and employment. L is the total available labor of the farmers in the two villages. The model assumes each farmer is endowed with the same acreage. The uppercase symbols denote aggregate-level quantities. |
$\begin{equation}
\pi_1^G = pb(l_f^G)^j(q_1^G + r_1)^ka_1^(1-j-k) + w(l_e^G) - [\gamma_s(1 - X_1) + \gamma_o(1 - X_2)] \\
\pi_2^G = pb(l_f^G)^j(q_2^G + r_2)^ka_2^(1-j-k) + w(l_e^G) - [\gamma_s(1 - X_2) + \gamma_o(1 - X_1)]
\end{equation}$ |
| The payoffs of group conformists in village 1 and 2. The cost of enforcement for a G increases with the frequencies of Os in their own village and the other village. |
$\begin{equation}
\pi_1^O = pb(l_f^O)^j(q_1^O + r_1)^ka_1^(1-j-k) + w(l_e^O) - \delta\bigg(1 - \sigma \frac{Q(L_m)}{Q(L_m*)}\bigg)q_1^O\frac{(X_1 + X_2)}{2}\\
\pi_2^O = pb(l_f^O)^j(q_2^O + r_2)^ka_2^(1-j-k) + w(l_e^O) - \delta\bigg(1 - \sigma \frac{Q(L_m)}{Q(L_m*)}\bigg)q_2^O\frac{(X_1 + X_2)}{2}
\end{equation}$ |
| The payoffs of opportunists in the two villages. The probability of opportunists being caught increases with the average frequency of Gs. The penalty varies by situation: it increases with the amount of water stolen, but decreases with water abundance in the system. |
$\begin{equation}
Q(L_m) = I(L_m)S
Q(L_m^*) = I(L_m^*)S
Abdundance = \frac{Q(L_m)}{Q(L_m^*)}
\end{equation}$ |
| Q(L-m*) is the optimal amount of water desired by the two villages. Q(L-m) is the total amount of water actually available in the system. |
$\begin{equation}
q_1^O \leq \frac{Q}{N_1} \\
q_1^G = min[\frac{Q}{N^G} , (Q - Q_1^O) / N_1^G]
\end{equation}$ |
| Due to their privileged access, Os in village 1 are less constrained by a higher bound on the amount of water they can steal. Gs in village 1 rely on their upstream position to extract water to bring their actual amount as close as possible to allocated amount. |
$\begin{equation}
q_2^O \leq \frac{Q - Q_1^O - Q_1^G)}{N_2}\\
q_2^G = \frac{[Q - Q_1^O - Q_1^G - Q_2^O]}{N_2^G}
\end{equation}$ |
| The amount of water that Os in village 2 can steal is constrained on the water available after consumption by Gs and Os in village 1. The second equation is the amount of water available for Gs in village 2. *Note that the model assumes that both Gs and Os in village 1 compete for water at the same time and they are symmetric in their capacity for competing for water. For the same reason, Os in village 2 face an upper bound on the water they obtain. |
#delta=2.9 / 1.0 R=30 / 0
par j=0.3,k=0.4
par p=1,b=1
par N1=50,N2=50
par ind_l=1
par ind_a1=1
par ind_a2=1
par Imax=1,psix=0.2, epsilonx=0.125, S=100
par w=0.2
par R=0
par gamma_s=0.05, gamma_o=0.1
par delta=1.4
par sigma=0.9
A1=(N1*ind_a1)
A2=(N2*ind_a2)
L=(N1+N2)*ind_l
ind_r1=R*(ind_a1/(A1+A2))
ind_r2=R*(ind_a2/(A1+A2))
psi=L*psix
epsilon=L*epsilonx
####################################################
#functions
infra(x)=if(x<=(psi-epsilon)) then (0) else (if(x>=(psi+epsilon)) then (Imax) else (((Imax)/(2*epsilon))*(x-psi+epsilon)))
f_income(lfx,qx,rx,ax)=p*b*((lfx)^(j))*((qx+rx)^(k))*((ax)^(1-j-k))
up_bound(x,up)=if(x>=up) then (up) else (x)
low_bound(x,low)=if(x<=low) then (low) else (x)
bound(x,low,up)=if(x<=low) then (low) else (if(x>=up) then (up) else (x))
step(x,ref)=if(x>ref) then (1) else (0)
####################################################
# first block
A=A1+A2
Qmax=Imax*S
CC=p*b*A^(1-j-k)
QR=(Qmax)/(2*epsilon)
####################################################
# Region 1 functions
Reg1Dec(Lf1tx)=if(Lf1tx<L) then (1) else (2)
Reg1Ymax(Dc,Y1max_ax,Y1max_bx)=if(Dc==1) then (Y1max_ax) else (Y1max_bx)
Reg1Lf(Dc,Lf1a_starx,Lf1b_starx)=if(Dc==1) then (Lf1a_starx) else (Lf1b_starx)
# Region 1
Lm1t=0
Lf1t=(CC*j*R^k/w)^(1/(1-j))
Lf1a_star=Lf1t
Y1max_a=CC*(Lf1t^j)*(R)^k+w*(L-Lf1t)
Lf1b_star=L
Y1max_b=CC*(Lf1b_star)^j*R^k
D1=Reg1Dec(Lf1t)
Lm_star_1=0
Lf_star_1=Reg1Lf(D1,Lf1a_star,Lf1b_star)
Y1max=Reg1Ymax(D1,Y1max_a,Y1max_b)
####################################################
# Region 3 functions
Reg3Dec(Lf3tx,Lm3tx,x)=if((Lf3tx+Lm3tx)<x) then (1) else (2)
Reg3Ymax(Dc,Y3max_ax,Y3max_bx)=if(Dc==1) then (Y3max_ax) else (Y3max_bx)
Reg3Lf(Dc,Lf3a_starx,Lf3b_starx)=if(Dc==1) then (Lf3a_starx) else (Lf3b_starx)
Reg3Lm(Dc,Lm3a_starx,Lm3b_starx)=if(Dc==1) then (Lm3a_starx) else (Lm3b_starx)
# Region 3
Lm3t=psi+epsilon
Lf3t=(CC*j*(Qmax+R)^k/w)^(1/(1-j))
Lf3a_star=Lf3t
Lm3a_star=Lm3t
Y3max_a=CC*Lf3a_star^j*(Qmax+R)^k+w*(L-Lf3a_star- Lm3a_star)
Lf3b_star=L-psi-epsilon
Lm3b_star=Lm3t
Y3max_b=CC*(Lf3b_star)^j*(Qmax+R)^k
temp=L
D3=Reg3Dec(Lf3t,Lm3t,temp)
Lm_star_3=Reg3Lm(D3,Lm3a_star,Lm3b_star)
Lf_star_3=Reg3Lf(D3,Lf3a_star,Lf3b_star)
Y3max=Reg3Ymax(D3,Y3max_a,Y3max_b)
####################################################
# Region 2 functions
Reg2Dec(L2sumx,Lm2tx,Lm2Dx,Lm2f_starx,Lf2f_starx)=if(L2sumx<=L) then (DecC(Lm2tx)) else (DecF(Lm2Dx,Lm2f_starx,Lf2f_starx))
DecB(Lm2tx)=if(Lm2tx<(psi-epsilon)) then (2) else (1)
DecC(Lm2tx)=if(Lm2tx>(psi+epsilon)) then (3) else (DecB(Lm2tx))
DecE(Lm2Dx)=if(Lm2Dx<(psi-epsilon)) then (5) else (4)
DecF(Lm2Dx,Lm2f_starx,Lf2f_starx)=if(Lm2Dx>(psi+epsilon)) then (DecG(Lm2f_starx,Lf2f_starx)) else (DecE(Lm2Dx))
DecG(Lm2f_starx,Lf2f_starx)=if((Lm2f_starx+Lf2f_starx)>=L) then (7) else (6)
Reg2Ymax(Dc,Y2max_ax,Y2max_bx,Y2max_cx,Y2max_dx,Y2max_ex,Y2max_fx,Y2max_gx)=if(Dc==1) then (Y2max_ax) else (if(Dc==2) then (Y2max_bx) else (if(Dc==3) then (Y2max_cx) else (if(Dc==4) then (Y2max_dx) else (if(Dc==5) then (Y2max_ex) else (if(Dc==6) then (Y2max_fx) else (Y2max_gx))))))
Reg2Lf(Dc,Lf2a_starx,Lf2b_starx,Lf2c_starx,Lf2d_starx,Lf2e_starx,Lf2f_starx,Lf2g_starx)=if(Dc==1) then (Lf2a_starx) else (if(Dc==2) then (Lf2b_starx) else (if(Dc==3) then (Lf2c_starx) else (if(Dc==4) then (Lf2d_starx) else (if(Dc==5) then (Lf2e_starx) else (if(Dc==6) then (Lf2f_starx) else (Lf2g_starx))))))
Reg2Lm(Dc,Lm2a_starx,Lm2b_starx,Lm2c_starx,Lm2d_starx,Lm2e_starx,Lm2f_starx,Lm2g_starx)=if(Dc==1) then (Lm2a_starx) else (if(Dc==2) then (Lm2b_starx) else (if(Dc==3) then (Lm2c_starx) else (if(Dc==4) then (Lm2d_starx) else (if(Dc==5) then (Lm2e_starx) else (if(Dc==6) then (Lm2f_starx) else (Lm2g_starx))))))
# basics
Lf2t=((j*CC/w)*(k*QR/j)^k)^(1/(1-j-k))
Lm2t=psi-epsilon-(R/QR)+(1/QR)*( (j*CC/w)*(QR*k/j)^(1-j) )^(1/(1-j-k))
L2sum=Lf2t+Lm2t
# Region 2
Lf2a_star=Lf2t
Lm2a_star=Lm2t
Y2max_a=(w/j)*(1-j)*Lf2a_star+w*(L-Lm2a_star)
Lf2b_star=(CC*j*R^k/w)^(1/(1-j))
Lm2b_star=psi-epsilon
Y2max_b=(w/j)*(1-j)*(Lf2b_star)+w*(L-Lm2b_star)
Lf2c_star=(CC*j*(Qmax+R)^k/w)^(1/(1-j))
Lm2c_star=psi+epsilon
Y2max_c=(w/j)*(1-j)*Lf2c_star + w*(L-Lm2c_star)
Lm2D=(k*L+j*(psi-epsilon)-j*R/QR)/(j+k)
Lf2D=L-Lm2D
Lm2d_star=Lm2D
Lf2d_star=Lf2D
Qd=R+QR*(Lm2d_star-psi+epsilon)
Y2max_d=CC*(Qd)^k*(Lf2d_star)^j
Lm2e_star=psi-epsilon
Lf2e_star=L-Lm2e_star
Y2max_e=CC*(L-Lm2e_star)^j*R^k
Lm2f_star=psi+epsilon
Lf2f_star=(CC*j*(Qmax+R)^k/w)^(1/(1-j))
Y2max_f=CC*(Lf2f_star)^j*(Qmax+R)^k+w*(L-Lm2f_star-Lf2f_star)
Lm2g_star=psi+epsilon
Lf2g_star=L-Lm2g_star
Y2max_g=CC*(L-psi-epsilon)^j*(Qmax+R)^k
D2=Reg2Dec(L2sum,Lm2t,Lm2D,Lm2f_star,Lf2f_star)
Lm_star_2=Reg2Lm(D2,Lm2a_star,Lm2b_star,Lm2c_star,Lm2d_star,Lm2e_star,Lm2f_star,Lm2g_star)
Lf_star_2=Reg2Lf(D2,Lf2a_star,Lf2b_star,Lf2c_star,Lf2d_star,Lf2e_star,Lf2f_star,Lf2g_star)
Y2max=Reg2Ymax(D2,Y2max_a,Y2max_b,Y2max_c,Y2max_d,Y2max_e,Y2max_f,Y2max_g)
####################################################
# Final decision functions
FinalDec(Y1maxx,Y2maxx,Y3maxx)=if(Y1maxx>Y2maxx) then (if(Y1maxx>Y3maxx) then (1) else (3)) else (2)
FinalYmax(Dc,Y1maxx,Y2maxx,Y3maxx)=if(Dc==1) then (Y1maxx) else (if(Dc==2) then (Y2maxx) else (Y3maxx))
FinalLm(Dc,Lm_star_1x,Lm_star_2x,Lm_star_3x)=if(Dc==1) then (Lm_star_1x) else (if(Dc==2) then (Lm_star_2x) else (Lm_star_3x))
FinalLf(Dc,Lf_star_1x,Lf_star_2x,Lf_star_3x)=if(Dc==1) then (Lf_star_1x) else (if(Dc==2) then (Lf_star_2x) else (Lf_star_3x))
# Final decision
DF=FinalDec(Y1max,Y2max,Y3max)
Ymax=FinalYmax(DF,Y1max,Y2max,Y3max)
Lm_star=FinalLm(DF,Lm_star_1,Lm_star_2,Lm_star_3)
Lf_star=FinalLm(DF,Lf_star_1,Lf_star_2,Lf_star_3)
Lo_star=L-Lm_star-Lf_star
aux Y_maxx=Ymax
aux Lm_starr=Lm_star
aux Lf_tarr=Lf_star
aux Lo_starr=Lo_star
###############################################
#======maintenance labor
lm_g1=min(Lm_star*(ind_a1/(A1+A2)), ind_l)
lm_g2=min(Lm_star*(ind_a2/(A1+A2)), ind_l)
lm_o1=0
lm_o2=0
#======Actual total maintenance labor
Lm=lm_g1*X1*N1+lm_g2*X2*N2
#======Actual water production
Q0=infra(Lm)*S
Qopt=infra(Lm_star)*S
Abund=Q0/Qopt
#======water allocation of upstream users
temp1=max( ((p*b*(ind_a1)^(1-j-k)*(j/w)^j*(k)^(1-j))/(delta*(1-sigma*Abund)*0.5*(X1+X2))^(1-j))^(1/(1-j-k))-ind_r1, 0 )
q_o1=min(temp1, Q0/N1)
Q1_steal=q_o1*(1-X1)*N1
q_g1=min(Q0*(1/(X1*N1+X2*N2)), (Q0-Q1_steal)*(1/(X1*N1)))
Q1_rule=q_g1*X1*N1
#======water left to downstream
Q1=Q1_rule+Q1_steal
Q2=max(Q0-Q1,0)
#======water allocation of downstream users
temp2=max(((p*b*(ind_a1)^(1-j-k)*(j/w)^j*(k)^(1-j))/(delta*(1-sigma*Abund)*0.5*(X1+X2))^(1-j))^(1/(1-j-k))-ind_r2, 0)
q_o2=min(temp2, Q2/N2)
Q2_steal=q_o2*(1-X2)*N2
q_g2=max(Q2-Q2_steal, 0)*(1/(X2*N2))
Q2_rule=q_g2*X2*N2
#======labor allocation of group conformist
lf_g1=min((p*b*j*((q_g1+ind_r1)^k)*(ind_a1^(1-j-k))/w)^(1/(1-j)), ind_l-lm_g1)
le_g1=max(ind_l-lf_g1-lm_g1, 0)
lf_g2=min((p*b*j*((q_g2+ind_r2)^k)*(ind_a2^(1-j-k))/w)^(1/(1-j)), ind_l-lm_g2)
le_g2=max(ind_l-lf_g2-lm_g2, 0)
#======labor allocation of opportunists
lf_o1=min((p*b*j*((q_o1+ind_r1)^k)*(ind_a1^(1-j-k))/w)^(1/(1-j)), ind_l-lm_o1)
le_o1=max(ind_l-lf_o1-lm_o1, 0)
lf_o2=min((p*b*j*((q_o2+ind_r2)^k)*(ind_a2^(1-j-k))/w)^(1/(1-j)), ind_l-lm_o2)
le_o2=max(ind_l-lf_o2-lm_o2, 0)
#======aggregate LF and LO
Lf=(lf_g1*X1+lf_o1*(1-X1))*N1+(lf_g2*X2+lf_o2*(1-X2))*N2
Le=(le_g1*X1+le_o1*(1-X1))*N1+(le_g2*X2+le_o2*(1-X2))*N2
#======upstream user income
f_g1=p*b*((lf_g1)^j)*((q_g1+ind_r1)^k)*(ind_a1^(1-j-k))
e_g1=w*le_g1
f_o1=p*b*((lf_o1)^j)*((q_o1+ind_r1)^k)*(ind_a1^(1-j-k))
e_o1=w*le_o1
#======upstream downstream income
f_g2=p*b*((lf_g2)^j)*((q_g2+ind_r2)^k)*(ind_a2^(1-j-k))
e_g2=w*le_g2
f_o2=p*b*((lf_o2)^j)*((q_o2+ind_r2)^k)*(ind_a2^(1-j-k))
e_o2=w*le_o2
#======net utility
pi_g1=max(f_g1+e_g1-(gamma_s*(1-X1)+gamma_o*(1-X2)), 0)
pi_o1=max(f_o1+e_o1-delta*q_o1*(1-sigma*Abund)*0.5*(X1+X2), 0)
pi_g2=max(f_g2+e_g2-(gamma_s*(1-X2)+gamma_o*(1-X1)), 0)
pi_o2=max(f_o2+e_o2-delta*q_o2*(1-sigma*Abund)*0.5*(X1+X2), 0)
#======right hand sides
pi_1=pi_g1*X1+pi_o1*(1-X1)
pi_2=pi_g2*X2+pi_o2*(1-X2)
dX1/dt=X1*(pi_g1-pi_1)
dX2/dt=X2*(pi_g2-pi_2)
#======aux
aux Q0x=Q0
aux v1_avg=(f_g1+e_g1)*X1+(f_o1+e_o1)*(1-X1)
aux v2_avg=(f_g2+e_g2)*X2+(f_o2+e_o2)*(1-X2)
aux v1_avgx=(f_g1)*X1+(f_o1)*(1-X1)
aux v2_avgx=(f_g2)*X2+(f_o2)*(1-X2)
aux q_g1x=q_g1
aux q_o1x=q_o1
aux q_g2x=q_g2
aux q_o2x=q_o2
aux q_o1_ax=temp1
aux q_o1_bx=q_o1
aux q_o2_ax=temp2
aux q_o2_bx=q_o2
aux Lfx=Lf
aux Lmx=Lm
aux Lex=Le
aux Lx=Lf+Lm+Le
aux lf_g1x=lf_g1
aux lm_g1x=lm_g1
aux le_g1x=le_g1
aux lf_o1x=lf_o1
aux lm_o1x=lm_o1
aux le_o1x=le_o1
aux lf_g2x=lf_g2
aux lm_g2x=lm_g2
aux le_g2x=le_g2
aux lf_o2x=lf_o2
aux lm_o2x=lm_o2
aux le_o2x=le_o2
aux f_g1x=f_g1
aux f_o1x=f_o1
aux f_g2x=f_g2
aux f_o2x=f_o2
aux e_g1x=e_g1
aux e_o1x=e_o1
aux e_g2x=e_g2
aux e_o2x=e_o2
aux pi_g1x=pi_g1
aux pi_o1x=pi_o1
aux pi_g2x=pi_g2
aux pi_o2x=pi_o2
aux pi_1x=pi_1
aux pi_2x=pi_2
yield=(1/p)*(f_g1*((X1)*N1)+f_o1*((1-X1)*N1)+f_g2*((X2)*N2)+f_o2*((1-X2)*N2))
ni=((f_g1+e_g1)*((X1)*N1)+(f_o1+e_o1)*((1-X1)*N1)+(f_g2+e_g2)*((X2)*N2)+(f_o2+e_o2)*((1-X2)*N2))
opt_yield=p*b*(L-Lm_star)^j*(Imax*S)^k*(A1+A2)^(1-j-k)
aux yieldx=yield
aux yieldxx=yield/opt_yield
aux nix=ni
s_f_1=(1/p)*(f_g1*((X1)*N1)+f_o1*((1-X1)*N1))/yield
s_f_2=(1/p)*(f_g2*((X2)*N2)+f_o2*((1-X2)*N2))/yield
gini_yield=abs(s_f_1-s_f_2)/(2*1^2*(s_f_1+s_f_2)/2)
aux gini_yieldx=gini_yield
s_ni_1=((f_g1+e_g1)*((X1)*N1)+(f_o1+e_o1)*((1-X1)*N1))/ni
s_ni_2=((f_g2+e_g2)*((X2)*N2)+(f_o2+e_o2)*((1-X2)*N2))/ni
gini_ni=abs(s_ni_1-s_ni_2)/(2*1^2*(s_ni_1+s_ni_2)/2)
aux gini_nix=gini_ni
aux Q0xx=Q0/100
#=============initial data
init X1=0.5, X2=0.5
@ bounds=10000
@ xp=X1, yp=X2
@ xlo=0, xhi=1, ylo=0, yhi=1
@ total=1000
@ dt=1
@ nmesh=90
@ XNC=1
@ YNC=7
@ MAXSTORE=1000000
done
#10=purple, 9=blue, 8=light blue, 7=green 6=light green, 5=gold, 4=orange, 3=darker orange, 2=red, 1=dark red, 0=blackVallury S, Arizona State University.
. 2015. Effect of infrastructure design on commons dilemmas in social−ecological system dynamics. Proceedings of the National Academy of Sciences. 112(43):13212.
